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Abstract 

Patriksson (Pat08l provided a then up-to-date survey on the continuous, separable, differentiable and convex 
resource allocation problem with a single resource constraint. Since the publication of that paper the interest in the 
problem has grown: several new applications have arisen where the problem at hand constitutes a subproblem, 
and several new algorithms have been developed for its efficient solution. This paper therefore serves three 
purposes. First, it provides an up-to-date extension of the survey of the literature of the field, complementing 
the survey in Patriksson IPat08: l with more then 20 books and articles. Second, it contributes improvements of 
some of these algorithms, in particular with an improvement of the pegging (that is, variable fixing) process in 
the relaxation algorithm, and an improved means to evaluate subsolutions. Third, it numerically evaluates several 
relaxation (primal) and breakpoint (dual) algorithms, incorporating a variety of pegging strategies, as well as a 
quasi-Newton method. Our conclusion is that our modification of the relaxation algorithm performs the best. At 
least for problem sizes up to 30 million variables the practical time complexity for the breakpoint and relaxation 
algorithms is linear. 


1 Introduction 

We consider the continuous, separable, differentiable and convex resource allocation problem with a single re¬ 
source constraint. The problem is formulated as follows: Let J := {1,2,..., n}. Let tpj : K. —> R and gj : R —>• K., 
j £ J, be convex and continuously differentiable. Moreover, let b £ K and — oo < lj < Uj < oo, j £ J. Consider 
the problem to 

minimize 0(x) := \ cf>j{xj), (la) 

jeJ 

subject to p(x) := X gj(xj) < b , (lb) 

jeJ 

lj < Xj < Uj, j £ J. (lc) 

We also consider the problem where the inequality constraint (fTbb is replaced by an equality, i.e.. 


minimize 

X 

<X x ) : = 

ieT 

(2a) 

subject to 

3( x ) : = X a i x i = 

j&J 

(2b) 


lj < Xj < Uj, j £ J, 

(2c) 
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where cij ^ 0, j £ J, and the sign is the same for all j £ J. Further, we assume that there exists an optimal 
solution to problems 0 and (0. For brevity, in the following discussions we define Xj := [lj,Uj], j £ J. 

Problems 0 and 0 arise in many areas, e.g., in search theory ( ||Koo99| ), economics ( 1Mar521 ). stratified 
sampling ([ BRS991 1. inventory systems ( 1MK93I 1. and queuing manufacturing networks (I BT89 I). Further, these 
problems occur as subproblems in algorithms that solve the integer resource allocation problem ( ]Mje83[ Sec¬ 
tion 4.7], BIK881 pp. 72-75], and lBS02bl ), multicommodity network flows ( !Sho851 Section 4.2]), and several 
others. Moreover, problems 0 and 0 can be used as subproblems when solving resource allocation problems 
with more than one resource constraint ( ]Mje83 [FZ8 31). and to solve extensions of problems 0 and 0 to a non- 
separable objective function </> ( |Mje83 l, DSV07 I); The books (Mje83|flK88i I.us )2 | describe several extensions, 
such as to minmax/maxmin objectives, multiple time periods, substitutable resources, network constraints, and 
integer decision variables. 

Many numerical studies of the problems 0 and 0 have been performed; for example, see I1BH81IINZ921 
[RJL921iKL981IKiw07HKiw08allKiw08bl . Our numerical study is however timely and well motivated, since except 
for those by Kiwiel I Kiw071 lKiw08al IKiw08b 1. where the quadratic knapsack problem is studied, none of the 
earlier approaches study large-scale versions of problem 0 or 0. There are also several algorithms (e.g., IINZ92I 
Section 1.4] and IISteQ 1II ) which are claimed to be promising, but have not been evaluated in a large-scale study. 
Only one earlier study ( IIKL981 ) evaluates the performance of algorithms for the problems 0 and 0 with respect 
to variations in the portions of the variables whose values are at a lower or upper bound at the optimal solution 
(see Section 16.21) . and this is done for modest size instances (n = 10 4 ) only. Further, no study has been done 
on the computational complexity for non-quadratic versions of the problems 0 or 0. Our numerical study 
also incorporates improvements of the relaxation algorithm, as presented in Sections 14.3.41 - 143751 and utilizes 
performance profiles ( IDM021 ). 

As a final note on the computational tests, we only consider problem instances where the dual variable cor¬ 
responding to the resource constraint (ITbl ). respectively (l2bl) . can be found in closed form; otherwise, we would 
need to implement a numerical method in some of the steps, e.g., a Newton method. We consider only customized 
algorithms for the problem at hand, since we presume that they perform better than more general algorithms under 
the above assumption. 

Patriksson jPat08l presents a survey of the history and applications of problems 0 and 0. Since its publi¬ 
cation several related articles have been published; the survey of l |Pat08 ;l is therefore complemented in Section[2] 
Section [3] presents a framework of breakpoint algorithms, resulting in three concrete representatives. Section Q] 
presents a framework of relaxation algorithms, and ultimately six concrete example methods. In Section 0 we 
describe a quasi-Newton method, due to Nielsen and Zenios 1NZ921 . for solving the problem 0. Section[6]de- 
scribes the numerical study. A test problem set is specified and the performance profile used for the evaluation is 
defined. In Section[7J we analyze the results from the numerical study. The structure is such that we first compare 
the relaxation algorithms, second the pegging process, and third the best performing algorithms among these two 
with the quasi-Newton method. Finally, we draw overall conclusions. 


2 Extension of the survey in [Pat08] 


We here extend the survey in IIPatQSl . using the same taxonomy, and sorted according to publication date. 

jMje83} £. M. Mjelde, Methods of the Allocation of Limited Resources, Section 4.7 

(Problem) <pj £ C 2 ; linear equality (dj = 1); lj = 0 
(Methodology) The ranking algorithm of |.I.(175 1 
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(Citations) Applications in capital budgeting (l !Han6811Shi771 ). cost-effectiveness problems ( tKir681[Pac70l 
|Mje78| ), health care ( lFet733 ). marketing ( 1LG7511 ). multiobjective optimization (|Geo673), portfolio 
selection ( lJdF75 il). production (the internal report leading to BBH79j ). reliability (|Bod69|), route¬ 
planning for ships or aircraft (I 1DBR661 ). search ( ICC58 1). ship loading ( ]Kyd69| ), and weapons selec¬ 
tion (Il )an67l ) 

(Notes) A monograph on resource allocation problems containing a comprehensive overview of the re¬ 
source allocation problem, including extensions to several resources, non-convex or non-differentiable 
objectives, integral decision variables, fractional programming formulations, etcetera. 

HSho85l N. Z. SHOR, Minimization Methods for Nondifferentiable Functions, Section 4.2 

(Problem) fj(xj) = ^(Xj — yj) 2 ', linear equality (aj = 1); lj = 0 

(Methodology) Pegging 

(Citations) ISI69L in which the motivating linear programming application is described 

(Notes) The problem arises within the framework of a right-hand side allocation algorithm for a large-scale 
linear program. 

HHZ05I z .-S. HUA AND B. ZHANG, Direct algorithm for separable continuous convex quadratic knapsack prob¬ 
lem (in Chinese) 

(Problem) <pj(Xj ) = ^ x 2 — VjXf, linear inequality (aj > 0); lj = 0 

(Methodology) Pegging 

(Citations) Algorithms for the problem ( I PK90; [MROOl |BS02bl |BS02al ) as well as for the case of integer 
variables 

(Notes) A numerical illustration (n = 6). 

ltDF06l! Y.-H. DAI AND R. Fletcher, New algorithms for singly linearly constrained quadratic programs sub¬ 
ject to lower and upper bounds 

(Problem) fj(xj) = x 2 — rjXj, qj > 0; gj convex in C 2 with g'(xj) > 0 

(Methodology) A combination of a bracketing algorithm on the Lagrangian dual derivative, and a secant 
algorithm for the Lagrangian dual problem 

(Citations) Algorithms for the problem ( I1HKL801 IBru84lICM871IPK901 ) 

(Notes) The problem arises as a subproblem in a gradient projection method for a general quadratic pro¬ 
gramming problem over a scaled simplex. 

ILS061 D. Li AND X. Sun, Nonlinear Integer Programming, Chapter 6: Nonlinear Knapsack Problems, Section 

6.1: Continuous-Relaxation-based Branch-and-Bound Methods 

(Problem) fj and gj increasing; tjj convex in C 2 with 9j>° 

(Methodology) Multiplier search 

(Citations) Multiplier search methods ( IIBS95I1 ). pegging methods ( iBS02bl[BS02a| ) 

(Notes) The problem arises as a subproblem in branch-and-bound methods for the integer programming 
version of the problem, such as for the quadratic knapsack problem, stratified sampling, manufacturing 
capacity planning, linearly constrained redundancy optimization in reliability networks, and linear cost 
minimization in reliability networks. 
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HDSV071 K. DAHIYA, S. K. SUNEJA and V. Verma, Convex programming with a single separable constraint 
and bounded variables 


(Problem) 4>j{xj) = ^fx 2 — TjXj, qj > 0; g 7 convex in C 2 with g'(xj) > 0; studies also the special case 
of a linear equality 

(Methodology) Iterative descent process using strictly convex quadratic separable approximations of a non- 
separable original objective / £ C 2 ; subproblems solved using the pegging algorithm of ISteOll 

(Citations) General references on convex programming over box constraints; 11 IKI.80 fDFL86l [PK901 for 
example algorithms for separable convex programming 

(Notes) Numerical QP (n = 6, g 7 quadratic) illustration; numerical comparison with an augmented La- 
grangian algorithm for a small problem (n = 2). 


tKiw07l K. C. KlWIEL, On linear-time algorithms for the continuous quadratic knapsack problem 

(Problem) <pj (xj ) = %fx 2 — TjXj, qj > 0; linear equality 

(Methodology) Breakpoint search algorithm applying median search of all breakpoints 

(Citations) Breakpoint search algorithms: | Bru84. , CM87 PK90 MSMJ03 I; sorting and searching meth¬ 
ods: lKnu98l[Kiw05l 

(Notes) Develops a general O(n) breakpoint algorithm; shows that the algorithms of | [PK90 MSMJ03 1 
may fail even on small examples; presents a modification of the breakpoint removal in the algorithm 
of I ICM87I . Numerical experiments (n £ [50 • 10 3 ,2 • 10 6 ]) for uncorrelated, weakly, and strongly 
correlated data; the new algorithm wins in CPU time over those in I Bru84 'l and ICM87I by 23%, and 
21%, respectively, on average. 


[ Kiw08al K. C. KlWIEL, Breakpoint searching algorithms for continuous quadratic knapsack problem 


(Problem) fj(xj) = x 2 — r 7 x 7 , qj > 0; linear equality 

(Methodology) A family of breakpoint search algorithms that include several choices of breakpoints for 
a median search and updates of quantities used for evaluating the piecewise linear implicit constraint 
function at the median point 


(Citations) Applications in resource allocation (I Bl 181 IBS97 , III 195 1). hierarchical production planning 
( IIBH81I ). network flows and transportation ( IIHKL80IISM90IIVen9 II INZ92IICH94II ), constrained ma¬ 
trix problems (I CDZ86 1), quadratic integer programming ([ BSS95 , BSS96 ,1II1951). Lagrangian relax¬ 
ation ( IHWC74I B. and quasi-Newton methods ( I1CM871 ); 0(n log n) sorting algorithms for the solu¬ 
tion of the Lagrangian dual problem (I HWC74i HKL801 ). 0(n) algorithms based on median search 
( 1Bru84l ICM87I IMdPJ89l |PK90l ICH94IIHH95IIMMP97I ) and approximate median search methods 
with 0(n) average-case performance ( HPK90I ); primal pegging algorithms with 0(n 2 ) worst-case per¬ 
formance (|Zip80l IBH8T1 IMic86l IVen91l IRJL92llBSS96l ) 


(Notes) Develops several variants of 0(n) breakpoint search algorithms, including some ideas earlier 
proposed in, e.g., I PK90. Cl 194 HH951 [MMP97 I; remarks that the more complex choices made in 
| MdP.I89.PK9(). 0194 HH95 i IMMP97 I also means that for some simple cases cycling may occur, and 
also provides convergent modifications for each of them. Numerical experiments (n £ [50 ■ 10 3 ,2 • 10 6 ]) 
for uncorrelated and weakly and strongly correlated data, and for both exact and inexact computations 
of the median; comparisons made with the 0(n) versions from I Brn84 CM87 I. reporting that a version 
(Algorithm 3.1) using exact medians is about 20% faster than the other ones; refers to an as yet un¬ 
available technical report from 2006 for more extensive tests and comparisons with pegging methods. 
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HKiw08bi K. C. KlWIEL, Variable fixing algorithms for continuous quadratic knapsack problem 


(Problem) &j(xj) = r 4fx’j — r ,Xj, qj > 0; linear equality 

(Methodology) Pegging 

(Citations) Applications: same references as in I IKiw08al : breakpoint search methods (1 HWC741IHKL801 
ICM87I lMdPJ89l IPK901ICH941IHH95IIMMP97IIMSMJ031 IKiw07l : pegging methods ( I1LG75I [BH791 
IBH8T1 ISho85l IMic86l IVen9 II IRJL92IIBSS96H 1 

(Notes) Develops a basic pegging algorithm and proposes several implementational choices for the solution 
of the reduced problem and the updates; shows that the algorithms in [ Mic861 IRJL92, , [BSS961 fail on 
a small example, and that there is a gap in the convergence analysis in |BH8l| (which also is closed) 
that affects algorithms whose analyses rest on that in I BII8I I (e.g., jVen9ll ); provide more efficient 
versions of several of these methods, including the introduction of incremental updates which reduce 
computations, and a more efficient stopping criterion. Numerical experiments (n £ [50 • 10 3 , 2 ■ 10 6 ]) 
for uncorrelated and weakly and strongly correlated data; comparisons made with the breakpoint search 
method of IIKiw08al which uses exact medians; on average the latter is 14% slower while at the same 
time it has a more stable run time; it is remarked that the advantage of pegging over breakpoint search 
has been reported also in IVen9HtRJL92 'l. 


itZHOSl B. ZHANG and Z. HUA, A unified method for a class of convex separable nonlinear knapsack problems 


(Problem) fij/g'j monotone and invertible, gl positive; consider Y^j=\gj( x j) <1 b where <] £ {<,=,>} 
(Methodology) Pegging algorithm using binary search on the value of fib/gb 


(Citations) Applications: resource allocation (11 .(175 Zip80 1BH8 1 , |Hoc94l ), the singly constrained multi¬ 
product newsvendor problem ( BHW63IIErl001lAM MM04 1), production/inventory problems ( 1BSSW94I 
IBSS95llBS02bll ). stratified sampling ( IIBSS95IIBRS99llBS02bll ). “core subproblem” f lUdF75IIAHKL80l 
IMT891 [RJL921IBSS95IIBSS961 ); algorithms: breakpoint methods ( llBS95l . itSteOTTl . IILG75I ) and re¬ 
laxation methods ( IIKL98II . I BS02bl . I BS02al ) 


(Notes) Claims (without a proof) that the complexity of the algorithm is 0(n), but the algorithm presented 
makes use of a mean-value method for the evaluation of the breakpoints which in the worst case results 
in 2n — 1 iterations. We also note that each iteration consist of ()(n) operations, whence the com¬ 
plexity is 0(n 2 ). Numerical experiments (n £ [10,10 4 ]) (fi quadratic, gj linear, <] chosen uniformly 
randomly). 


BDWW12II A. De WAEGENAERE AND J. L. Wielhouwer, A breakpoint search approach for convex resource 

allocation problems with bounded variables 

(Problem) fij and gj convex; strictly monotone and gj ([(/>' /<?j-] _1 ) is either strictly increasing or strictly 
decreasing for all j 

(Methodology) Breakpoint search algorithm using a refined pegging method (5-sets pegging); generalizes 
the quadratic breakpoint algorithm in IIPK90II and its extension in | Kiw08a] such that it applies for the 
problem setting in (BS02bl 

(Citations) Generalizes the quadratic breakpoint algorithm in I PK90I I and its extension in [ Kiw08a ! such 
that it is valid for fj and gj as in lBS02bl . Applications in resource allocation IIPK901INZ921 IBS95] 
|yW95llBSS96llBS02bllDW'()2llCLZ091IBSSVQ6llPat()8l 

(Notes) Discuss their algorithms’ advantages compared to other articles presented in IPat081 . 
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0KW12I G. KlM and C.-H. Wu, A pegging algorithm for separable continuous nonlinear knapsack problems 
with box constraints 

(Problem) </>' invertible: linear equality (aj > 0) 

(Methodology) Pegging; improves the pegging algorithm in I BH811 by allowing the primal box constraints 
to be checked implicitly, using bounds on their Lagrange multipliers 

(Citations) Applications to portfolio selection IIMar521 . multicommodity network flows 1AHKL80 1, trans¬ 
portation BOK84B . production planning I ITamSOI 

(Notes) Compares their methodology with the method in I BH81 1 on random test problems. On randomized 
quadratic continuous knapsack problems (n £ [5 • 10 3 ,2 • 10 6 ]) their algorithms wins by 8-10%, 
except for the smallest problems where the method in I BH81 1 wins by 12.5%. Two other types of 
test problems are also investigated, wherein the quadratic continuous knapsack problem arises as a 
subproblem: quadratic network flows, and portfolio optimization. In the former case, the algorithm 
(based on conjugate gradients) is taken from l AraOO I. Here, the pegging algorithm proposed wins over 
that in llBH8ll by 10-48%, with n £ [200, 5 ■ 10 3 ]. In the portfolio optimization algorithm, which 
is based on a progressive hedging method described in f AraOOfl , the quadratic continuous knapsack 
problem arises as a subproblem. Here, the range of n is not completely disclosed; however, the speed¬ 
up over the method in IBH811 is reported to be 21-25%. 

OLuslH H. LUSS, Equitable Resource Allocation: Models, Algorithms and Applications , Chapter 2: Nonlinear 
Resource Allocation 

(Problem) <f>j strictly convex; gj(xj) = Xj 

(Methodology) Pegging 

(Citations) |i Koo53l as an origin; I KL98I for computational examples; jPat()8 l as a survey 

(Notes) This book extends the resource allocation problem (discussed only in Chapter 2) in several ways, 
including equitable optimization through the use of minmax/maxmin objectives, multiple time periods, 
substitutable resources, network constraints, and integer decision variables. 

lfZCl2il B. ZHANG and B. Chen, Heuristic and exact solution method for convex nonlinear knapsack problem 
(Problem) fj strictly convex; gj(xj ) = Xj 

(Methodology) The problem at hand is a subproblem in a branch-and-bound procedure for the solution of 
an integer-restricted version of the problem 

(Citations) Applications to the newsvendor problem ( BErlOO , AMMM04B ), resource allocation ( IBH81I , 
IHoc941 ). production ( ||BSS95 1). and stratified sampling (II BRS991 ); efficient methods for the continu¬ 
ous relaxation (l BS951[KL981[Ste0l [ ; ZH083 ); heuristics for the integer program based on rounding of 
the solution to the continuous relaxation ( IBS951 lHZL06 iP; algorithms for the integer program based 
on the solution of continuous problems and branch-and-bound (l lBS95l[BS02bl l 
(Notes) Utilizing the algorithm from IIZH08 I1 to solve the continuous relaxations (and rounding to pro¬ 
duce feasible solutions), the authors develop a branch-and-bound algorithm. It is compared with the 
methods from jBS95l as well as with branch-and-bound algorithms utilizing a variety of tree search 
principles, on instances with quadratic objectives, according to problem generation principles from 
BBS02b> (n £ {10,15,500,1000,2000}). 

BBGRS13I L. Bayon, J. M. Grau, M. M. Ruiz, AND R M. SUAREZ, An exact algorithm for the continuous 
quadratic knapsack problem via infimal convolution 
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(Problem) <p 3 strictly convex quadratic; linear equality ( a 3 = 1) 

(Methodology) Sorting of breakpoints 

(Citations) Previous algorithms for the problem (( MT931ICH941 ) 

(Notes) Two numerical applications: (1) an economic dispatch problem (n = 5), and an academic example 
(n € [200,10 4 ]); in the latter example the results are favourably compared with the Matlab solver 
QUADPROG. 

[ DHH13I T. A. DAVIS, W. W. Hager, and J. T. HunGERFORD, The separable convex quadratic knapsasck 
problem 

(Problem) 4>j(xj) = x? — r 3 Xj ; q 3 > 0; linear equality 

(Methodology) Breakpoint search utilizing a heap data structure, based on an initial multiplier estimate 
using a secant Newton method 

(Citations) Applications HIHKL80IICM871ISM901INZ921IDF06I ); multiplier algorithms ( IIHKL80I lBru84l 
ICM87llPK90llMSM.T031 1: pegging methods r iBH81IISho85llMic86ll\^n91IIRJL92llBSS96llKiw08bl , l: 
quasi-Newton methods (I DF06 . CMSI4 I) 

(Notes) Numerical experiments (n = 6.25 • 10 6 ) on random test problems, examining (a) the best number of 
initial (secant) Newton iterations, and (b) the performance against a primal pegging algorithm. Overall, 
breakpoint search is favourable on its own given a good initial mulitiplier estimate; otherwise, 3 or 4 
iterations of the secant Newton method is a good initialization procedure. 

IFGI5I A. FRANGIONI AND E. GORGONE, A library for continuous convex separable quadratic knapsack prob¬ 
lems 

(Problem) = %fxj — rjXj\ q 3 > 0; gj(xj) = Xj 

(Methodology) Compares a breakpoint algorithm with CPLEX and concludes that the breakpoint algorithm 
outperforms CPLEX 

(Citations) Applications in resource allocation and algorithms ( IPat08l ) 

(Notes) Presents an open source library for the continuous convex separable quadratic knapsack problem 
and concludes that the library can be useful for further studies of the problem at hand. 

nzccsi3n t. Zhu, W. CHEN, J. Chen, and W. Sun, Direct algorithm for continuous separable knapsack 
problem (in Chinese) 

(Problem) <f>j(xj) = — r :) x 3 ; q 3 > 0; g 3 (xj) = a 3 Xj 

(Methodology) Pegging 

(Citations) Algorithms for the more general integer version of the problem; algorithms by Bretthauer and 
She tty ( IIBS95I lBS02bl lBS02al l: other specialized algorithms for the problem ( HBH8 1 1 TPK901 [RJL92I 
ISteOll IKiw08bl 1 

(Notes) A detailed numerical example (n = 8); favourable comparisons with a Matlab solver (n £ 
{50,100,200}). 

HZhal3l L. ZHANG, A Newton-type algorithm for solving problems of search theory 

(Problem) cpj(xj) = —a,j( 1 — exp(-CjXj)), a 3 ,Cj > 0; X is a scaled unit simplex 


7 






























(Methodology) The KKT conditions, as defined in ([3]l. are relaxed into a system of non-smooth equations 
through the utilization of the Fischer-Burmeister ( IIFis92l ) smoothing function; a Newton-like algo¬ 
rithm is then employed for these equations for a sequence of values of the smoothing parameter. The 
algorithm is shown to asymptotically and superlinearly converge to the unique optimal solution 

(Citations) Survey ( lPat08l ); methodologies ( lSte041 1: applications ( IIKoo99l ) 

(Notes) Three sets of numerical experiments. Main experiment (n £ [10 2 ,10 4 ]) on randomized problems; 
compares with the pegging algorithm in iSte04l . noting that the proposed algorithm is faster. The 
proposed methodology is however terminated based on a nonzero value of the Fischer-Burmeister 
smoothing function, whence the final solution need not be feasible or optimal. Second set of experi¬ 
ments on a problem taken from IWG691 (n = 4), showing no comparisons. Third experiment on data 
from the Bureau of Water Conservatory (n = 5), showing no comparisons. 

l iCMS14j R. COMINETTI, W. F. MASCARENHAS, AND P. J. S. Silva, A Newton’s method for the continuous 

quadratic knapsack problem 

(Problem) <f>j strictly convex quadratic; linear equality 

(Methodology) An approximative Newton method for the Lagrangian dual problem utilizing a secant glob¬ 
alization and variable fixing 

(Citations) Applications to resource allocation ( BBH81I III 195 , BS97I ), multicommodity network flows 
(| HKL80 i : NZ92ilSM90 3), Lagrangian relaxation using subgradient methods ( IHWC74IB . quasi-Newton 
updates with bounds ( ICM871 ), semismooth Newton methods ( 1FM04ID ; related methods ( 1DF06I 
|Kiw08b |), where the latter method is shown to be equivalent to the proposed one when there are only 
lower bounds or lower bounds 

(Notes) Numerical experiments (n £ [50 • 10 3 , 2 • 10 6 ]) comparing the proposed method to a secant method 
( BDF06I ). breakpoint search ( 1Kiw08bl ). and median search ( jKiw08al on problems with uncorrelated, 
weakly correlated, and correlated data. The proposed Newton method is overall better—about 30% 
better on the larger instances. A further test is made on the classification problems described in 1DF06II 
arising in the training in support vector machines; as the Hessian is non-diagonal, a projected gradient 
method is used, leading to subproblems of the type considered. Here, Newton’s method is superior. 

IWR141 S.E. WRIGHT and J. J. Rohal, Solving the continuous nonlinear resource allocation problem with an 

interior point method 

(Problem) fj and gj convex, and in C 2 on an open set containing [lj,Uj], Further, <pj is decreasing on 
[lj,Uj\ and gj is increasing on [lj,Uj], and Ylj^j 9j(h) < b < Y^j^j 9j( u j)- Test instances include 
resource renewal \0j(xj ) := CjXj(e^ 1 ^ Xi — 1) for Xj > 0, 0j(xj) := —Xj for Xj < 0, and gj linear], 
weighted p-norm over a ball \<pj(xj) := Cj(xj — yj) p and gj(xj) := \xj\ r with p,r £ {2, 2.5,3,4}], 
sums of powers [0j(x 3 ) := Cj\x 3 — yj\ Pj , and Pj{xj) := |£jj ri ], convex quartic over a simplex [0j a 
fourth-power polynomial, and g 3 (xj) := Xj], and log-exponential [fj{xj) := exp(aijXj+dij)], 

gj linear] 

(Methodology) A damped feasible interior-point Newton method for the solution of the KKT conditions 

(Citations) Survey (j Pat08l ); application to resource renewal (| MR00 | ); methodologies ( |Bru84[ |PK90l 
IRJL921IMR0Q1 IKiwOSall 1 

(Notes) The algorithm is introduced for problems where subsolutions are not available in closed form. 
Shows that the linear system defining the Newton search direction is solvable in 0{n ) time. Numerical 
experiments (n £ [10 2 ,10 6 ]) conclude that the interior point method wins over breakpoint search, often 
by an order of magnitude. 
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3 Breakpoint algorithms 

Algorithms based on the Lagrangian relaxation of the explicit constraint ( ITbt have an older history than the re¬ 
laxation algorithms. This is probably due to the fact that the relaxation algorithm quite strongly rests on the 
Karush-Kuhn-Tucker (KKT) conditions, which did not become widely available until the end of the 1940s and 
early 1950s with the work of F. John HJoh48fl . W. Karush |Kar39l . and H. W. Kuhn and A. W. Tucker KEE 
Lagrangian based algorithms have been present much longer and the famous “Lagrange multiplier method” for 
equality constrained optimization is classic in the calculus curriculum. Indeed, Lagrange multiplier techniques for 
our problem <(T]) date back at least as far as to flBec521 ; see RPat08B for a survey of the history of the problem. 

Considering problem <[TJ and introducing the Lagrange multiplier /t for constraint ( I I h| i. we obtain the follow¬ 


ing conditions for the optimality of x* in <jT]) : 

F*>0, ff(x*)<6, n*(g(x*)-b) = 0, (3a) 

x* G X j: j € J, (3b) 

and 

x*=l jt if c/)j(xj) > j G J, (3c) 

x* = u jt if </>'(£*)<j G J, (3d) 

lj < X* < Uj, if fijiXj) = ~n*g'j{x*), j G J. (3e) 


For a fixed optimal value /i* of the Lagrange multiplier the conditions (I3cl>-(l3el> are the optimality conditions for 
the minimization over x G TTjLi Xj of the Lagrangian function defined on n; =1 Xi x R + , 

n 

L(x, n) := -bn + ^2{<t>j(xj) + ngj(xj)}- 

i=i 

Given n > 0 its minimization over x G YYj=i Xj separates into n problems, yielding the Lagrangian dual function 

n 

q{n) := — bfj, + y minimum {<pj(xj) + ngjixj)}- (4) 

* ^ Xj€Xj 

3 =1 

By introducing additional properties of the problem, we can ensure that the function q is not only concave but 
finite on M + and moreover differentiable there. Suppose, for example, that for each j, (j)j (■) + ng :) (•) is weakly 
coercive on Xj for every g > 0 [that is, that either Xj is bounded or that for every /r > 0, <pj (xj) + n9j i x j) tends 
to infinity whenever Xj tends to ±oo], and that 0 :) is strictly convex on Xj. In this case the derivative q' exists on 
R + and equals 

'/( f ) = + n9j( x j(v)), 

where x(/i) is the unique minimum of the Lagrange function L(-, /i) over rij=i Xj. Thanks to this simple form of 
the dual derivative, the maximum /i* of q over R + can be characterized by the complementarity conditions (l3al) . 
and the conditions (0 are the primal-dual optimality conditions for the pair of primal-dual convex programs. 

If we assume that /i* ^ 0, we search for n* > 0 such that q'{n*) = 0 [° r > m other words, g(x(/r*)) = 6], that 
is, we need to solve a special equation in the unknown entity /i, where the function q' is implicitly defined, but is 
known to be decreasing since q is concave. This equation can of course be solved through the use of any general 
such procedure [for example, bisection search takes two initial values Ji and /i with q'(]T) < 0 and q'{n) > 0, then 
iteratively cancels part of the initial interval given the sign of q' at its midpoint (Jl + /r)/2], but the structure of q' 
makes specialized algorithms possible to utilize. 
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From the above optimality conditions for the Lagrangian minimization problem, we obtain that 

f lj, if ft > ■= 

x Av) = l Uj, if H<ltf ■- -<j/ j {u j )/g , j {u j ), j G J. ( 5 ) 

if <t>'j {Xj ) + fig'jixj) = °, 

In a rudimentary algorithm we order these indices (or, breakpoints) g l j and in an increasing (for example) order 
into {hi, ..., //y}, where N < 2n due to the possible presence of ties. Finding /i* then amounts to finding an 
index j* such that that q'(g.j*) > 0 and q '< 0; then we know that g* G g 3 *+i) and q'(g*) = 0. 
Hence from equation (0, we know for all j if x* = lj , x* = Uj or lj < x* < Uj. Now by fixing all variables x* 
that equal the corresponding lower or upper bound we can ignore the bound constraint ( fTB and find an analytical 
solution of the problem. 

Two decisions thus need to be made: how to find the index j*, and how to perform the interpolation. Starting 
with the former, the easiest means is to run through the indices in ascending or descending order to find the index 
where q 1 changes sign. If we have access to the indices j + and for which q , (/r J +) > 0 while q\g 0 ~) < 0, then 
we can choose the midpoint index, check the corresponding sign of q', and reduce the index set accordingly. Given 
the sorted list, we can also find this index in some randomized fashion. 

As remarked above, algorithms such as bisection search can be implemented without the use of the break¬ 
points, and therefore without the use of sorting, as long as an initial interval can somehow be found; also general 
methods for solving the equation q'(g) = 0, such as the secant method or regula falsi, can be used even without an 
initial interval; notice however that q (f_ C 1 , whence a pure Newton method is not guaranteed to be well-defined. 

While the sorting operation used in the ranking and bisection search methods takes 0(n log n) time, it is 
possible to lower the complexity by choosing the trial index based on the median index, which is found without 
the use of sorting; the complexity of the algorithm is then reduced to 0(n). It is not clear, however, that the latter 
must always be more efficient, since the “O” definition calls for n to be “large enough”. 

We also remark that in the case when the problem ([T]) arises as a subproblem in an iterative method, as the 
method converges the data describing the problem will tend to stabilize. This fact motivates the use of reoptimiza¬ 
tion of the problems, which most obviously can be done by using the previous value of the Lagrange multiplier as 
a starting point and/or utilizing the previous ordering of the breakpoints; in the latter case, the 0(n log n) sorting 
complexity will eventually drop dramatically. 

In Section l3Tl we consider the breakpoint algorithm for the equality problem 0. In Section lT2l we describe 
three pegging methods and in Sections 13.3113.51 we apply these pegging methods to the breakpoint algorithm. 
Finally, in Section [T6l we briefly discuss the convergence and time complexity of the breakpoint algorithm. 

3.1 Equality constraints 

We now consider problem 0 where the inequality of the primal constraint is replaced by an equality. For the 
problem to be convex, the resource constraint (l2bl > has to be affine, i.e., g(x) := (l j x j ~ b. Beside the 

resource constraint, the Lagrangian and the optimality conditions will take the same form as for problem 0 but 
with one important difference; // is unrestricted in sign, whence the condition < f3al > is replaced by “g(x*) = 6”. 

3.2 The pegging process 

The origin of the pegging process is found in the relaxation algorithm; see, e.g., | BH8Tj |. The purpose of pegging 
is to predict properties of the primal variables in the optimal solution from an arbitrary dual value. In Sections 
13.2.1113.2.21 and !3.2.3l we show how to predict if an optimal primal variable value equals its lower or upper bound, 
or is strictly within any of the bounds. 
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3.2.1 2-sets pegging 

If we can determine if a variable equals its lower bound at the optimal solution then we add its variable index to a 
set L and reduce the original problem. Similarly, if we know that a variable equals its upper bound at the optimal 
solution then we add the variable index to the set U. Using the sets L and U when solving problem © or © will 
be referred to as 2-sets pegging. 

Assume that we have a lower limit // and an upper limit Ji on the optimal dual value, that is, // < /r < Ji. From 
© we can define the sets L(g) := {j £ J | /£ > —<f>j{lj)/g'j{lj)} and U(JL) := {j £ J \ Jl < —<ftj{uj)/g'j{uj)}. 
Let J k ■= J\ {L(g) U U{-p)} and let b k := b - Ej £ l( h ) 9Ah) ~ E je u(ji) 9j{uj)- Hence we can define a 
subproblem of problem dTJ as follows: 


minimize <^(x) '■= E (6a) 

j&J k 

subject to p(x) := E 9 j(xj) < b k , (6b) 

j£J k 

lj < Xj < Uj, j G J k . (6c) 

Similarly we can define a subproblem of problem © as follows: 

minimize 0(x) := E (7 a) 

j£j k 

subject to p(x) := ^ UjXj = b k , (7b) 

jeJ k 

lj < Xj < Uj, j € J k ■ (7c) 


Consider problem ©. Assuming that p* > 0, the constraint ( ITbt has to be fulfilled with equality. For any given 
dual variable ///'' we can determine the primal solution x /r (jj k ) of problem © from condition ©. We know from 
Section [3] that all optimality conditions © except the resource constraint l ITbl ) are satisfied. Moreover, we know 
that the resource constraint has to be fulfilled with equality, in order for the solution to be optimal. Substituting x fc 
into the resource constraint will be referred to as explicit evaluation; this leaves us with three cases, namely 


E 9j( x j) = hk i 

(8a) 

j<EJ k 


E 9A x j ) <b k ,or 

(8b) 

j£j k 


E aA xk ) > >>'"■ 

(8c) 

j£J k 



If (f8ab is fulfilled for x A: then all optimality conditions are met and x* = x ,:: . Consider next the case ( l8bl >; clearly 
x fc is not optimal but we know that x* is such that EjeJ^ 9j( x j) > 9j( x j) si nce Ejej fc 9j( x j) = b k - The 

function <j 3 is convex and differentiable and can increase in one interval and decrease in another (consider, e.g., 
gj{xj ) = x'j), which implies that no predictions can be made of the size of x* relative to that of x k . Hence, we 
need gj to be monotone. For problem ©, Bretthauer and Shetty | BS02bl Section 2] consider four cases equivalent 
to the following: 

Case 1: For all j e J, g 3 is decreasing and n{xj) := —<j>'j{xj)/g'j(xj) is increasing in Xj. 

Case 2: For all j £ J, gj is increasing and jj{x :] ) is decreasing in x 3 . 
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Case 3: For all j £ J, g 3 is decreasing and g(xj) is decreasing in x 3 . 

Case 4: For all j £ J, g 3 is increasing and gix 3 ) is increasing in x :j . 

If Case 3 or 4 holds it is possible to find a closed form of the optimal solution to the problem 0, see I BS02h 
Proposition 10], Considering problem 0, Cases 3 and 4 cannot occur, since the resource constraint is affine, i.e., 
g(xj) := —(f>'j(xj)/aj. Hence, only Cases 1 and 2 are of interest here. 

Note that if g(xj ) is increasing in xj then x 3 (/./.) is nondecreasing and vice versa. This is an essential property 
for the validity of the pegging process. We can indeed state the following proposition (a similar one can be stated 
for problem (0, but without the assumption g* > 0): 

Proposition 1 (pegging for Cases 1 and 2). Consider problem 0 and assume that g* > 0. 

(i) If Case 1 holds, and if < [8bt holds, then x* = lj for all j £ L(g k ). 

(ii) If Case 1 holds, and if (l8cl > holds, then x* = u 3 for all j £ U{g k ). 

(iii) If Case 2 holds, and if (l8bl > holds, then x* = Uj for all j £ U{g k ). 

(iv) If Case 2 holds, and if (l8cl > holds, then x*j = lj for all j £ L(g k ). 

Proof. A proof of (i) is given; the proofs for (ii), (iii), and (iv) are analogous. From (l8bl > we have that 

9j{xj{p k )) < b k . 

jeJ k 

In Case 2, for all j gj is increasing and Xj(g) is nonincreaing, which implies that gj(xj(g)) is nonincreasing in g 
for all j. Hence, we have that g k = JI > g* which implies that x k < x* for all j since Xj(g) is nonincreasing in 
g for all j. Hence, for j £ U ( g k ) we know that x k+1 = Uj = x*. i.e., we can peg j £ U ( g k ). □ 

3.2.2 3-sets pegging 

As in the 2-sets pegging principle of Section [3.2.1l we determine if a variable takes the value of the lower or upper 
bound at the optimal solution. Additionally for the 3-sets pegging we determine if a variable belongs to the open 
interval between the lower and upper bound, i.e., if x* £ ( lj,Uj ). Assume that we know that g.* £ (/i“, g l j); then it 
follows from 0 that x* £ (l :l . u :] ) and there is no need to check if x k equals the lower or upper bound which will 
reduce future calculations. 3-sets pegging is used for the quadratic knapsack problem (f24l> in IIKiw08al Section 3]. 
The method described in [ Kiw08al Section 3] can be generalized according to the following proposition: 

Proposition 2 (relax primal variables from lower and upper bounds). Assume that Case 1 or 2 in Section 13.2.11 
holds and that for some values of g and JI, g < g.* < Ji holds. If g,,JI £ [g. 1 ^ /r“] holds for some j £ J k then 

lj < Xj < Uj. 

Proof. Assume that Case 2 holds, i.e., g = —fb/g'j is decreasing. (The proof for Case 1 is analogous.) Assume 
that /i, JI £ [/i“, g l j holds for some j £ J k . Since — </>' / gl is decreasing we have that g,JI £ [gfj , g/j\ implies that 
g* £ (g'j, g l j ), and from 0 it then follows that lj < x* < Uj. □ 

3.2.3 5-sets pegging 

As in the 3-sets pegging principle in Section l3.2.2l we determine if a variable takes the value of the lower or upper 
bound or if the variable strictly belongs to the interval between the lower and upper bound. For 5-sets pegging 
we also determine if a variable is larger than the lower bound or smaller than the upper bound. 5-sets pegging for 
problem 0 is used in 0DWW12I. generalizing a method from [fKiw08al. Assuming that we know that X* < Uj, 
there is no need to check if x k equals the upper bound; this might reduce future calculations. The proof of the 
following proposition follows from the monotonicity of g 3 and x 3 (ji) (see IDWW12 1 for a proof). 
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Proposition 3 (relax primal variables from lower or upper bound). Assume that Case 1 or 2 in Section [3. 2.1 I holds. 
Let j £ J k . If pt* < pt < n l j, then x*j > lj and if pt* > pt > fij then x* < Uj. 


□ 


3.3 Algorithm: Median search of Breakpoints with 2-sets pegging (MB2) 

Consider Case 1 or 2 in Section l3.2.1l for problem 0. Let median(-) denote a function which provides the median 
of a finite vector; let pt m be the median breakpoint, and define the total use of the resource due to variables that 
equal the lower and upper bounds as ft := EjefAV. < fc } 9j (h)’ and Pu : = T,je{N\^>^ m } 9j K0> respectively. 

In the spirit of PBS02al Section 3.1] and | Kiw071 Algorithm 3.1], we present the following algorithm: 

Initialization: 

Set N := J,k := 1, and b k := b. 

Compute breakpoints pt ; := (pt* ) jeN , /*“ := (/iV) j6 n as in <l5j- and let p i k := (—oo, (ft*) 1 , (pt“) T , oo) T . 
Step 0 (check if /i = 0 is optimal): 

If E je jv 9j ( x j (0)) — b, for x 3 (pt) is determined from (0, then pt* = 0, x* = Xj (0) for j € N. Stop. 

Iterative algorithm: 

Step 1 (stopping test): 

If pt fc = 0 then find x* and pt* from problem 0 relaxed from lower and upper bounds. 

Otherwise, let pt m := median(pt k ). 

Step 2 (compute explicit reference): 

Determine <5 := J2je{N\^<^ m <^ 1 .} 9j{ x j{Vm)) + Pu + Pi, where x 3 (fi) is determined from ©. 

If <5 > b k , then go to Step 3.1. 

If <5 < b k , then go to Step 3.2. 

Otherwise (5 = b k ) let pt* := pt m , find x* from 0, and stop. 

Step 3.1 (update and fix lower bounds): 

For all j G N: If pt m > /jJ 3 then let N := N \ {j} and x* := lj. 

Let n k+1 := {n k j ) je { N \ tlm <^}, b k+1 := b k - fa, and k := k + 1. 

Go to Step 1. 

Step 3.2 (update and fix upper bounds): 

For all j £ N: If pt m < pt“ then let N := N \ {j} and x* := u 3 . 

Set pt fc+1 := {9j)je{N\n m >n^}, b k+1 := b k - fi u , and k := k + 1. 

Go to Step 1. 

For problem 0, the algorithm is similar except Step 0 vanishes. 

3.4 Algorithm: Median search of Breakpoints with 3-sets pegging (MB3) 

If we can determine if the value of a variable x 3 is strictly within the bounds for all pt such that // < // < pt, then we 
don’t have to check if x 3 violates the bounds when we determine x 3 from 0 in future iterations (see Proposition 
0. This might save us some operations. Define a set of indices for the lower and upper limit being within the 
interval of the lower and upper breakpoint for a variable: M := {j £ N \ fi,~p £ [fj l 3 , pt“] }. From Proposition [2] 
we have that if j £ M then lj < x* < Uj. Hence if j £ M we do not have to check if x 3 violates the bounds in 
future iterations. But we should have in mind that to determine if j £ M requires some extra operations. If we let 
the initial values of the limits be pt = —oo and pi = oo, then ji. ~p (j [pt*,pt“] and we note that there is no need to 
check if j £ M, as long as pt = — oo or pt = oo. Hence, to avoid unnecessary operations we start by checking if 
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j £ M when // > — oo or JJ < oo; this is in contrast to the algorithms in IKiw08al Section 3] and IDWW12 1. 

We hnally define the contribution to the resource constraint from the variables including M : 7 (fi) := 
Y,eM Note that the value of 7 depends on the value of the parameters in <pj, gj and fi. If we can 

separate the parameters from the multiplier fi, i.e., if 7 is additively and/or multiplicatively separable [that is, 
7 (fi, ■) = 7 i(/tt) 72 ( - ) + 73 (’)] then the values of 72 and 73 can be calculated successively so that no calculations are 
done more than once. Consider, for example, the negative entropy function, <f>j{xj) := Xj log(a 'j/a,j — 1) and the 
resource constraint function gj(xj) := x :j . Then, 7 ^, a) := YjeM 9j( x j (m)) = YjeM a j e ~ M = 72 (a) 7 i(/x), 
where 71 (/r) = e _M and 72 (a) = YjeM a v i- e -’ we update 72 in Steps 3.1 and 3.2 such that if //, /T £ [//“, g}-\ 
then 72 := 72 + a :) . For the quadratic knapsack problem, this approach is applied in IIKiw08al Section 3], We next 
present an algorithm that applies the usage of M; Steps 0 and 1 are similar to the algorithm MB2 in Section [PI 

Initialization: 

Set N := J, k := 1, b k := b, 7 := 0, M := 0, fi := — 00 , and Jl := 00 . 

Compute breakpoints fi 1 := (fij)j e N, /<“ := (fi’j)jeN us in (0, and let fi k := 

Iterative algorithm: 

Step 2 (compute explicit reference): 

Determine S := Yje{N\^<^ m <^} 9j( x j(Hrn)) + 7(/7n) +/3 U + Pi, 
where Xj (fi) is determined from 0 . 

If S > b k , then go to Step 3.1. 

If S < b k , then go to Step 3.2. 

Otherwise (<5 = b k ) set fi* := fi m \ find optimal x* from 0, and stop. 

Step 3.1 (update and fix lower bounds): 

Let fi — fi m . 

For j £ N: If fj, > fi l j then let N := N\{j} and x* := lj. 

If fi,Ji € [fij, fi l j\ then let N := N \ {j} and M := M U {j}; update 71 and 72 . 

Let fi k+1 := (fij) je{NUM \^ m <^}, bk+1 ■= bk - A. and k := k + 1. 

Update 7 and go to Step 1. 

Step 3.2 (update and fix upper bounds): 

Let fi . — fim- 

For j £ N: If ~p < fij then let N := N \ {j} and x* := Uj. 

If h,T l £ ■ l l 'j] then let N := N \ {j} and M := M U {j}; update 71 and 72 . 

Let fi k+1 := (k-j) je{ NuM\^ m >^}, bk+1 ■= b k - Pu, and k := k + 1. 

Go to Step 1. 

Remark 1. If 7 is not separable then the updates of 71 and/or 72 in Steps 3.1 and 3.2 vanish. 

3.5 Algorithm: Median search of Breakpoints with 5-set pegging (MB5) 

As in the algorithm MB2 in Section 13.31 we peg the variables whose values equal the lower or upper bounds, 
and as in MB3 we determine if a variable is strictly within the bounds. Further, we determine if a variable is 
smaller than the upper bound or larger than the lower bound as in Section 13.2.31 Define a set L _ of indices 
where the optimal solution is known to be strictly smaller than the upper bound in the optimal solution, i.e., 
L_ := {j £ N | x* < Uj}, and similarly define a set U + of indices where the variable is known to be larger than 
the lower bound in the optimal solution, i.e., U+ := {j £ N \ x* > lj}. Define, respectively, the total use of the 
resource due to variables whose values equal the lower and upper bounds as /3i := Yje{NuL \n l <v } //j ((j ) anc * 



14 










f} u := Ylje{Nuu+\iJ. u >n } ( u j )• We present an algorithm that applies 5-sets pegging where Steps 0 and 1 are 

similar to the algorithm MB2 in Section l3T3l 

Initialization: 

Set N := J, k := 1, b k := b, 7 := 0, and M = L_ = U+ := 0, /i := — 00 , and Jl := 00 . 

Compute breakpoints fj, 1 := (lJ l ,j)jeN, n u '■= (iij)j^N as in @ 1 , and let /x k := 

Iterative algorithm: 

Step 2 (compute explicit reference): 

Determine 6 := Eje{Nui/+UL^<^ m <^} 9j{%j(Mm)) + l(Mm) + Pu + Pi, 
where Xj(fi) is determined from (|5). 

If 5 > b k , then go to Step 3.1. 

If 8 < b k , then go to Step 3.2. 

Otherwise (5 = b k ) set fj* := and find optimal x* from (0, and stop. 

Step 3.1 (update and fix lower bounds): 

Let [A . = fj>m 

For j G N: If n > n l j then let N := N \ {j} and x* := lj. 

For j 6 L_: If ^ > fi l j then let := L_ \ { j } and x* := lj. 

For j G N: If /i > /x“ then let N := N \ {j} and L_ := L_ U {j}. 

For j e U + : If n G [/r“, fil] and JI > ji'- then let U + := U + \ {j} and M := M U {j}; update 7 l5 72 . 
Let n k+1 := {^j) je{NUM uL.uu+\fi m <^}’ bk+1 ■= bk ~ ^ and k := k + 1. 

Go to Step 1. 

Step 3.2 (update and fix upper bounds): 

Let .— fiiji. 

For j G N: If Jl < /r“ then let N := N \ {j} and x* := Uj. 

For j G U + : If JI < /r“ then let U + := U+ \ {j} and x* := Uj. 

For j G N: If Jl < jjl then let N := N\{j} and U+ := U+ U {j}. 

For j G L-\ If/T 7 [/z“, fj,j\ and /i < ^ then let L- := L- \ {j} and M := M U {j}; update 71 , 72 . 
Let /r fc+1 := (Vj) je{NUM L) L _uu + \nm>^}’ bk+1 := ^ ^ and k ~k+ 1. 

Go to Step 1. 

Note that when we determine x 3 (i.i rn ) for j G U + we do not need to check if x 3 = l 3 , and for j G L_ we do need 
to check if x 3 = Uj. This will in some cases save us some operations. Further when we determine if x 3 G M we 
only need to check if this holds for j G U + or j G L_ depending on if we update the lower or upper bound of the 
dual variable. This differs from the algorithm in I DWW12 1 since that algorithm finds M from U+ U L_. Note that 
we make use of information from earlier iterations when updating <5. When updating 6 in IDWW121 information 
from earlier iteration is negligible hence some operations are repeated. 

Remark 2. Consider iteration k. In Step 3.1 we only need to check if fi m G [//"', /Jj) and JJ > f.i'- for j G U+ since 
we know from earlier iterations that if j G N then Jl ^ fi l 3 . Similar holds for Step 3.2. 

3.6 Convergence and time complexity of breakpoint algorithms 

Similar to jKiw08al and | I)WW12 | it is possible to show that the breakpoint algorithms converge to the optimal 
solution. Consider the time complexity for algorithm MB2, MB3 and MB5. Assuming that the median function 
median(-) is linear, we have operations in each of the Steps 0 through 3.2 for some constants Ci for i G 
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{0; 1; 2; 3.1; 3.2} corresponding to Step 0 through 3.2 respectively. Since we use a median search function the 
number of iterations will terminate in log 2 (n) iterations (in the worst case). Hence the time complexity of the 
algorithm is 0[n log 2 (n)). However in HBru84l a proof for a time complexity of 0(n) is given for the breakpoint 
algorithm solving the quadratic knapsack problem (l24l >. 

4 Relaxation algorithms 

In a relaxation algorithm for the problem 0 we iteratively solve the problem relaxed from constraints (flcl) . i.e., 
we solve the following problem: 

minimize </>(x) := E (9a) 

jeJ k 

subject to g(x) := ^ 9j( x j) < b. (9b) 

jejfc 


From problem (0 we obtain a solution x. Together with x we also obtain an estimate p of the multiplier value //* 
from the optimality condition. Then we adjust the solution x for constraints (ITcl) by determining Xj from 

! lj > if Xj < lj, 

Uj , if Xj>Uj, (10) 

Xj , if lj < Xj < Uj . 

At the beginning of the algorithm we must determine whether constraint (flbl) is satisfied with equality at an 
optimal solution. From the optimality condition (l3al) we have that if the inequality constraint ( I I hi is satisfied then 
fj* = 0. Hence, for p = 0 we find Xj by solving 4>j(xj) + pgj(xj) = 0; if Yhjej 9j( x j) — b, where the value 
of Xj is determined from (ITOt . then we have found the optimal solution to problem 0. Otherwise, we know that 
fi* > 0 and that the inequality constraint (ITbl) can be regarded as an equality. Hence, we solve the problem 0, 
obtaining a solution x. Let 


L(x) := {j £ J k | Xj < lj }, U(x) := {j € J k \ Xj > Uj } 

denote the sets of variables that are either out of bounds at x or equal a lower, respectively an upper, bound. 

In order to simplify the remaining discussion, we consider Case 2 in Section 15.2.11 i.e., n(xj) to be mono- 
tonically decreasing and <j :j is increasing; Case 1 is treated analogously. Calculate the total deficit and excess at x, 
respectively, as 

V := E^ife) - 9j(xj)), A := ^2(gj(xj) - 9j{uj)). (11) 

jeL jeu 

Now, if A > V then we set x* = Uj for j e ?7(x); if A < V we set x* = lj for j £ L(x); otherwise A = V and 
we have found the optimal solution. If A 7 ^ V, then we reduce the problem by removing the fixed variables, and 
adjusting the right-hand side of the constraint (ITbt to reflect the variables fixed. If any free variables are left, we 
re-solve problem 0 and repeat the procedure, otherwise we have obtained an optimal solution. 

The rationale behind this procedure is quite simple and natural: Suppose that A > V holds. We have that 
p = —<j)'j{xj)/g'j{x) for j G J k \{IU U}. Let s £ U(k) and i £ J k \ t/(x). Since the functions — 0 ' / g' ] are 
decreasing, it follows that 

<t>' 8 {u>s) > <j>' s {Xs) = . = fijiXj) > <t>'j{Ui) 

g'sM - 9s(x a ) 9i(xi) - 9j(ui)' 
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Denote by b + the right-hand side in the following iteration given that A > V holds: b + := b — f2j<=u{x) 9j( u j)- 
Also let (x' ; p') denote a pair of relaxed optimal primal-dual solutions in the following iteration. We must have 
that p! < p, since 


55 ft to) = b- Y 9j to) <b- Y ft to) =b+ = 

j&J k \u(x) jeu(x) 


j£U(x) 


55 

j£j k \U(x) 


hence, for at least one j £ J\U (x) we have that x' > Xj, and therefore. 


M = 


<t>j(x'j) , 


< — 


9'Ax'-) g'Ax) 


= P 


follows. This derivation was first described by Bitran and Hax BBH811 . 

Since in each iteration at least one variable is pegged to an optimal value, the algorithm is clearly finite. The 
most serious disadvantage of the algorithm may be the requirement that the problem without the variable bounds 
present must have an optimal solution. The computational efficiency of this method is also determined by whether 
or not it is possible to provide an explicit formula for each Xj in terms of the multiplier. 


4.1 Explicit/Implicit evaluation 

For the breakpoint algorithm we determine x from (0 for an arbitrary multiplier //; then we explicitly evaluate the 
optimality of x by substituting it into the resource constraint ( ITbt . The explicit evaluation leaves us with one out 
of three possible scenarios: (f 8 ab— (f 8 cl> . For the relaxation algorithm the traditional method to evaluate a solution to 
the problem © is to calculate the total deficit and excess at x (see Section [Hi. Also this evaluation leaves us with 
3 possible scenarios, namely 


A = V, 

( 12 a) 

A > V, 

( 12 b) 

A < V. 

( 12 c) 


Evaluating the optimality from (1 1 2ab - (l 1 2cl ) will be referred to as an implicit evaluation. For the relaxation 
algorithm we next show that the implicit and explicit evaluations are equivalent. Propositions [4] and 0 below 
state the relations between explicit and implicit evaluation. The proof of Proposition [His similar to the proof of 
Proposition^ 

Proposition 4 (relation between explicit and implicit evaluation for g :j monotonically increasing). If for all j £ J 
gj is monotonically increasing, then the explicit evaluation (f8al)—(f8cl> is equivalent to the implicit evaluation (U2al >- 
(112cl) . i.e.. (112al) ([Sat . (112bl) -£=> (l8bl) . and (112cl) <£=>■ (l8cl >. 

Proof. For all j £ J , let x :J be the solution to ©. Let V and A be defined as in (ITTI) . We have from (fTOt that 

X! % to) = 51 * to) + Y to fe) - ft to) + ft to)} 

j&J jGJ\{UUL} j£L 

+ Y ift to) ~ ft to) + fi'ito)} 

jeu 

= Y ft to) + Y {ft (to - 9j to)} + Y ft to) 

j£j\{UUL} j£L j£L 

- Y {ft to) - 9j to)) + Y ft to) 

jeu jeu 

= 55 ft to) + v - A 

= b + V- A, 
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where the last equality follows from the fact that x is the solution to the relaxed problem. Hence, x must satisfy the 
resource constraint. We know that A, V > 0 since g :j is increasing. Hence, if A = V then (l8al) holds, if A > V 
then ( l8bl > holds, and if A < V then (l8cl > holds. □ 

Proposition 5 (relation between explicit and implicit evaluation for gj monotonically decreasing). If g :/ is mono¬ 
tonic ally decreasing for all j € J, then the explicit evaluation (f8ab— (f8ct is equivalent to the implicit evaluation 
(l8at—(l8c1>. i.e.. (1 1 2ab <*=> (Hal l. (1 1 2cb <*=>- (f8bt . and (1 12bb -£=>■ (Hel l. 

4.2 Primal/Dual evaluation of boundaries 

In the breakpoint algorithms in Section [3] we use the dual variable to determine if Xj equals a bound or if it lies in 
between the bounds; see (0. In relaxation algorithms the traditional way to solve the problem 0 is to determine 
the primal optimal solution Xj of the relaxed problem 0 and then simply check if Xj is within the lower and upper 
bounds; see (ITOt . An alternative is to find the optimal dual variable ft of the relaxed problem 0 and then (similar 
to the breakpoint algorithm in Section 0 determine the primal variables from 0. Of course this requires us to 
determine the breakpoints. However, for the variables that violate the bounds (ITcl) we don’t have to determine x 7 
from the relation 4>j(xj) + gg 3 (xj) = 0. Hence, instead of evaluating L(x) and (7(x) from the primal variable as 
in Section[4] we can evaluate L(p) and U(p.) from the dual variable as in Section [3.2.1l 

Define <!>., (xj) '■= 4>'j ( x j ) / 9j ( x j ) an d assume that there exists an inverse to <1^. It follows from the optimality 
conditions 0 that = —p,, anc * ®j{ u j) ~ ( u j ) / 9'j( u .i )• Hence we have that: 

■= {j e J | A > -4>'j{h)/9'j{h )} <=> £(*) ■= {j e J | % < l 0 }, 

and 

U{p) ~ {j G J | A < -<t>j(v>j)/g'j{uj) } U(k) := {jeJ\ xj > u 3 }. 

4.3 Implementation choices for the relaxation algorithm 

Considering the performance of a relaxation algorithm two decisions need to be made: should the relaxed solution 
of the problem be evaluated implicitly from V and A or explicitly from 0 g :1 (see Section 14.11 )7 Should the 
algorithm solve the primal or dual relaxed problem (see Section l4~2l i? An overview of the relaxation algorithm 
and its possible realizations is shown in Figure Q] The leftmost path in the figure is the classic primal relaxation 
algorithm of Bitran and Hax I1BH81I . the rightmost path in the figure is implemented in l lSteOl i. and the algorithm 
that applies the right path in Step 1 and the left path in Step 2 is implemented in ||KW12| . Beside these two paths 
no other paths have been explored. Our intention is to evaluate the theoretically and practically best performing 
paths. Since no earlier studies have applied 3- or 5-sets pegging for the relaxation algorithm, our intention is to 
apply these two more sophisticated pegging methods. 

In Section 14.3.1 1 we present an algorithm corresponding to the leftmost path in Figure Q} in Section [03] we 
present an algorithm which utilizes the rightmost path in Step 1 of Figure [I] and then changes to the left path in 
Step 2 of the figure; in Section l4.3.3l an algorithm corresponding to the rightmost path of Figure|T]is presented; and 
in Section l4.3.4l an algorithm that utilizes the rightmost path in Step 1 of FigureQ]and then utilizes the theoretically 
best path in Step 2 is presented. All algorithms in Sections [4. 3. 1144. 3T4l utilize 2-sets pegging. In Section l4.3.5l we 
describe how 3- and 5-sets pegging can be utilized in these. 

4.3.1 Algorithm: Primal determination with Implicit evaluation of the Relaxed problem (PIR2) 

We first assume that Case 1 in Section [3.2.1l holds for problem 0. Define parameters to calculate the total deficit 

and excess, as a\ := 9j( x j ). a k _ := E jeL * 9j( x j), /3 k := gj{lj) and := Yjeu” 9j(uj)- 
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Step 1 


Step 2 


Step 3 



FIGURE 1 : The diagram shows different possibilities to consider when constructing a relaxation algorithm. 


Hence we have that A fc = on. — dt and V A ' = /3_ - a_. In the next iteration when we reduce the resource b k we 
hence don’t have to re-calculate the part of the pegged variables which define or /3— (see Steps 3.1 and 3.2). 
We present an algorithm for problem (Q} which is similar to the algorithm in Section 2 of I BS02bl : 

Initialization: Set k := 1, J k := J, and b k := b. 

Step 0 (check if /i = 0 is optimal): 

Let n = 0 and find the solution x to the relaxed problem ([9]», i.e., solve 4>'j(xj) + fig'j{xj), j £ J k . 
If E, e Jk 9 j{xj) < b, for Xj determined from (ITOt . then /i* = 0, x*j = Xj for j £ J k , and stop. 
Iterative algorithm: 

Step 1 (solve relaxed primal problem): 

For j £ J k , find x k by solving the relaxed problem (|9). 

Set L := 0 and U := 0. 

Step 2 (implicit evaluation): 

Determine (7(x fc ) and L(x k ) while computing A k := a* — ,3+ and V fc := 3- — &-■ 

If A k > V fc , then go to Step 3.1. 

If A k < V k , then go to Step 3.2. 

If A k = V fc , then set x* := lj for j £ L(x fc ), x* := Uj for j £ C/(x fc ), 
x* := x k for j £ J k \ {L(x. k ) U C/(x fc )}, and stop. 

Step 3.1 (peg lower bounds): 

Set x* := lj for j £ L, b k+1 := b k - /3 k and J k+1 := J k \ L. 

If J k+1 := 0 then stop, else set k := k + 1 and go to Step 1. 

Step 3.2 (peg upper bounds): 

Set x* := Uj for j £ U, b k+1 := b k - 3\ and J k+1 := J k \ U. 

If J k+1 := 0 then stop else set k := k + 1 and go to Step 1. 

We need to clarify some of the steps of the algorithm. In Step 1, we find x fc+1 from, or partly from, x fc . Assume, 
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for example, that <j)j (xj ) = Xj log(xj/aj — 1) and Qj (xj) = Xj ; then, x k+1 = ajb k+1 / ^ J ej fc + 1 a j =b k+1 /(uj — 
Yhj£K{it k ) a j)’ where u> = Yl jeJ k a J anc * K :=U if the upper bound was pegged at iteration k and K := L if the 
lower bound was pegged at iteration k. If \K\ < \J k+1 | then this will save us some operations. A similar update 
of x for the quadratic knapsack problem is performed in | ;RJL92l Section 3], I BSS96 I and I KiwOSh Section 5.1]. 

As in IIKiw08bl Algorithm 3.1], our algorithm will stop if A k = V /:: , while the algorithm in | BS02h Section 
2] stops only if L(x k ) U U ( x k ) = 0. Moreover, in Steps 3.1 and 3.2, we peg the variables that violate the bounds 
and calculate b k explicitly, while in IBS02bl Section 2] the index j is added to the set of violated bounds (L or U) 
and b k is calculated as b - - Eject 9j( u i)- 

According to Proposition [T] if Case 2 in Section 13.2.1 1 holds. Step 2 in the algorithm is modified as follows 
(an analogous modification can be defined for the relaxation algorithms in Sections [4.3.2l through l4.3.5l >: 

Step 2’ (implicit evaluation): 

Determine (7(x fc ) and L(x fc ) while computing A k = a+ — /3+ and V fc = /3^ — a* . 

If A k > V fc , then go to Step 3.2. 

If A fc < V fc , then go to Step 3.1. 

If A fc = V fc , then set x* := lj for j £ L(Sc k ), x* := Uj for j £ U(ic k ), 
x* := x k for j £ J k \ (i(x fe ) U t/(x fc )}, and stop. 

Remark 3. For the equality problem (0, // is unrestricted; hence the algorithm for problem 0 will be similar to 
the above, except we ignore Step 0. 

Remark 4. In Step 1 we have to calculate Xj \-J k \ times. In Step 2 we need to find Xj which needs at most 2\ J k \ 
comparisons. We also have to calculate and A fc , which implies 2| L U U\ operations. 

4.3.2 Algorithm: Dual determination with Implicit evaluation of the Relaxed problem (DIR2) 

Instead of evaluating the primal variables in Step 1 as in the algorithm PIR2 in Section 14.3.11 we can evaluate 
the dual variable /i. Note that if we evaluate the dual variable, then we have to determine the breakpoints in the 
initialization. Our modification of the algorithm in Section l4.3.1l takes the following form (Steps 0 and 3 are same 
as in PIR2 and are therefore not repeated): 

Initialization: 

Set k := 1, J k := J, and b k := b. 

Calculate breakpoints // := (jJ l j)j(zjk, and fi u := as in 0. 

Step 1 (solve relaxed dual problem): 

Find the optimal dual variable fl k of the relaxed problem 0. 

Step 2 (implicit evaluation): 

Determine L(jj, k ) and U(fi k ) while computing A k := a\ — and V fe := /3^ — a* . 

If A k > V fc , then go to Step 3.1. 

If A k < V fc , then go to Step 3.2. 

If A fc = V fc , then set x* := lj for j £ L(fi k ), x* := Uj for j £ U(p, k ) 
x* := x k for j £ J \ {L{fi k ) U U(fi k )}, and stop. 

Remark 5. In Step 1 we need to find the optimal dual solution fj, k to the relaxed problem. In Step 2 we need at 
most 21 J k \ comparisons but we only need to calculate Xj(fi) for j £ {L U U}. For the evaluation we need to 
calculate V fe and A k , which implies 2|L U U\ operations. 
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4.3.3 Algorithm: Dual determination with Explicit evaluation of the Relaxed problem (DER2) 


As in the algorithm DIR2 in Section 14.3.21 we evaluate the dual variable but instead of evaluating the relaxed 
solution implicitly we do it explicitly. Define Pi := 9jih) and Pu ~ IZjeu 9j( u j)- The algorithm follows 

(Steps 0 and 3 are same as in PIR2 in Section [4.3. 1 l and are therefore not repeated): 

Initialization: 

Set k := 1, J k := J, and b k := b. 

Calculate breakpoints p i 1 := fi u := as in (0. 

Step 1 (solve relaxed dual problem): 

Find the optimal dual variable fr of the relaxed problem (0. 

Step 2 (explicit evaluation): 

Determine U (/t fe ) and L(fi k ) and calculate 5(fi k ) := YjjeJ k \{U{p,*)uL{p,*)} 9j(xj(P k )) + Pi + Pu¬ 
li 5(p k ) > b k , then go to Step 3.1. 

If 5(fi k ) < b k , then go to Step 3.2. 

If S(p k ) = b k , then set x* := lj for j £ L(/t fe ), x* := Uj for j G U(p k ) 
x* := x k for j G J \ {L((i k ) U U(p k )}, and stop. 

The algorithm uses the principle of the algorithm in ESteOll Algorithm 1], 

Remark 6. In Step 1 we need to find the value of jx k from the relaxed problem. In Step 2 we need at most 2| J k \ 
comparisons and we need to calculate x k for j £ J k \{L U U}. For the evaluation we need to calculate S(fi k ), 
which implies |J fc | operations. 

4.3.4 Algorithm: Dual determination modification with blended evaluation of the Relaxed problem (DBR2) 

Consider the implicit evaluation in Section 14.3.21 We calculate V fc and A fc from V A: := 9j(^j) ~ 

E je r/(A fc )ffi(^ fc ) and Afc : = E ieC /(A fe ) 9j(Xj) - 9j(uj), which implies 2P\L k U U k \ operations, 

where P is an integer associated with the number of operations it takes to calculate gj(xj). Moreover, we have to 
determine Xj for j £ {L k U U k }, which implies Q\L k U U k \ operations, where Q is an integer associated with the 
number of operations it takes to determine Xj{(i). 

Now consider the explicit evaluation: We have to calculate 5{jl k ) := Ej<= j fe \{[/(A fe )uL(A fe )} 9j( x j{p k )) + 
Ejei(A fe ) 9j(h) + Ej6t/(A fe ) which implies P\J k \ operations. Moreover, we have to determine Xj for 

j £ J k \ {L k U U k }, which implies Q\J k \ {L k U U k }\ operations. 

Hence, if (P + Q)\ J k \ < (2 P + 2Q)\U(^ k ) U L(/r fc )| or, equivalently, | J k \ < 2|C7(/z fc ) U L(/i fc )|, then using 
the explicit evaluation of the relaxed solution x in Step 2 would require less operations, and it would be more 
successful to use the algorithm in Section 14.3.31 If however | J k | > 2\U(fi k ) U L(/i fc )|, then there will be less 
operations if we use the algorithm in Section l4.3.2l So, we propose a new algorithm that utilizes the cardinalities 
of the sets J k , U ( p k ) and L(p k ): from the cardinalities we make the decision whether to use an explicit or implicit 
evaluation in Step 2. We consider the following modification of the algorithm PIR2 in Section l4.3.1l 
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Initialization: 

Set k := 1, J k := J, and b k := b. 

Calculate breakpoints fi l := := {k L 'j)j^J k as i n ©• 

Step 1: (solve relaxed dual problem): 

Find the optimal dual variable ji k of the relaxed problem ©. 

Step 1.1: (implicit or explicit evaluation?): 

Determine U(jl k ) and L{(l k ). 

If | J k \ < 21 U{(i k ) U L(p, k ) then use explicit evaluation (continue with algorithm DER2 in 
Section l4.3.3b . otherwise use implicit evaluation (continue with algorithm DIR2 in Section l4.3.2b . 


4.3.5 3- and 5-sets pegging for relaxation algorithms 

From the proof of convergence of the relaxation algorithm (see Section l4~4l > we have that the algorithm improves 
the lower or the upper bound for the dual variable in each iteration. Hence, similar to the breakpoint algorithm (see 
Sections [3.41 and [375] ) it is possible to apply 3- and 5-sets pegging for the dual relaxation algorithms DIR2, DER2 
and DBR2. Similar to the relaxation algorithm we will denote a relaxation algorithm that uses 3-sets-pegging with 
a suffix ”3”, e.g., DBR3, and one that uses 5-sets-pegging with a suffix ”5”, e.g., DBR5. The implementation of 
3- and 5-sets pegging is similar for the three dual relaxation algorithms; therefore only DIR3 and DIR5 are given 
(Step 0 is similar as for PIR in Section l4.3.1l >: 

DIR3: 

Initialization: 

Set k := 1, J k := J, M k := 0, b k := b, /z = — oo, and Jl = oo. 

Calculate breakpoints := {ji l j)j e jk, \i u := as in (0. 

Step 1 (solve relaxed dual problem): 

Find the optimal dual variable ji k of the relaxed problem 0. 

Step 2 (implicit evaluation): 

Determine L(fi k ) and U(p, k ) from J k and compute A fc = a\ — 0+ and V fc = — a^_. 

If A k > V fc , then go to Step 3.1. 

If A k < V fc , then go to Step 3.2. 

If A k = V fc , then set x* = lj for j e L(fx k ), x* = Uj for j GU ( (x k ) 
x* = Xj(p, k ) for j G {J k U M k } \ {L(p, k ) U (7(/t fc )}, and stop. 

Step 3.1 (peg lower bounds): 

Update lower bound, [i := fi k , and resource b k+1 := b k — (3^. 

Set x* := lj for j G L k , J k+1 := J k \ L k and M k+1 := M k . 

For j G J k+1 : If /£,/! G then J k+1 := J k+1 \ {j}, M k+1 := M k+1 U {j}. 

If J k+1 = 0 then find optimal solution and stop, else set k :=/;■+ 1 and go to Step 1. 

Step 3.2 (peg upper bounds): 

Update upper bound, jl := jl k , and resource b k+1 := b k — 0i_. 

Set x* := Uj for j G U k , J k+1 := J k \ U k and M k+1 := M k . 

For j G J k+1 : If ft,Ji G [/.then J k+1 := J k+1 \ {j}, M k+1 := M k+1 U {j}. 

If J k+1 = 0 then find optimal solution and stop, else set k := k + 1 and go to Step 1. 

In Steps 3.1 and 3.2, if J k+1 = 0 then we can find the optimal solution from M since we know that for all j G M 
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it holds that lj < x :j < Uj, i.e., we don’t have to consider the constraint die] ). Further, in Step 3.1 (and analogusly 
for Step 3.2) when searching for j £ J k+1 such that /j. Jl £ [/x“, //*■], we only have to consider j £ J k+1 \ U k , 
since we know that if j £ U k then ji k < //". Note that for the case where ji k = fij the corresponding index j 
will continue to belong to J k+1 . Finally we note that we don’t have to check if /i, Jl £ [/i“, /rj] if /x = — oo and/or 
Jl = oo. 

DIR5: 

Initialization: 

Set k := 1, J k := J, M k = U\ = L k _ := 0, b k := b, /£ = -oo, and Jl = oo. 

Calculate breakpoints /j/ := := {l^])jeJ k as i n •13- 

Step 1 (solve relaxed dual problem): 

Find the optimal dual variable Ji k of the relaxed problem <[9]). 

Step 2 (implicit evaluation): 

Determine L(/i fc ) and U(Ji k ) from J k U L k _ U U k , compute A k = a\ — /3+ and V fe = /3^ — a*. 

If A k > V k , then go to Step 3.1. 

If A k < V k , then go to Step 3.2. 

If A k = V k \ then set x* = lj for j £ L(fx k ), x* = Uj for j £ U (/ l k ) 

x* = X j{Jj, k ) for j £ {J k U M k UL^U U'j} \ {L(jl k ) U U{Jx k )}, and stop. 

Step 3.1 (peg lower bounds): 

Update lower bound, /j, := Ji k , and resource b k+1 := b k — j3_. 

Set x* := ^ for j £ L k , J k+1 := J k \ L k and L k _ +1 := L k _ \ L k . 

For j £ {J k+1 \ U k }\ If Hm < ji) then J k+1 := J k+1 \ {j} and L k _ +1 := L k + X U {;/}. 

For j £ {U% \ U k }: If/x < fi] then Ul := Ul \ {j} and M k := M k U {j}. 

Set M k+1 := M k and [/+ +1 := U k . 

If J k+1 U L k i +1 U Ul +1 = 0 then find optimal solution and stop, 
else set k := k + 1 and go to Step 1. 

Step 3.2 (peg upper bounds): 

Update upper bound, jl := ji k , and resource b k+1 := b k — 

Set x* := Uj for j £ U k , J k+1 := J k \ U k and U k+1 := U k \U k . 

For j £ {J k+1 \ L k _}: If /j rn < //' then J k+1 := J k+1 \ {y} and [7^ +1 := U^ +1 U {j}. 

For j £ {L k _ \ L k }: If fj, < then L k _ := L k _ \ {j} and M k := M k U {j}. 

Set M k+l ~ M k and L k _ +1 := L k _. 

If J k+1 U L* +1 U U k+1 = 0 then find optimal solution and stop, 
else set k := k + 1 and go to Step 1. 


Remark 7. Note that for DBR2 in Section 14.3.41 we determine whether we should continue with DIR or DER 
from | J k \ < 2\U(jl k ) U L{jl k )\. This condition needs to be modified for the 3- and the 5-sets algorithms DBR3 
and DBR5. This is done by substituting the condition mentioned by \J k U M k \ < 2\U k (jl k ) U L k (jl k )\ and 
| J k U M k U L k _ U U^l < 2\U k (jl k ) U L k (jl k )\ for 3- and 5-sets pegging, respectively. 

4.4 Optimality of the relaxation algorithms 

For the inequality problem (Q}- optimality and validation for the pegging process for the primal relaxation algorithm 
PIR2 was established by Bretthauer and Shetty (BS02b, Propositions 1-9], Let k* be the iteration where the 
algorithm terminates. Then Bretthauer and Shetty state that the primal relaxation algorithm generates the following 
solution for problem ( 0 : 
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Bretthauer and Shetty establish that this solution fulfills all KKT conditions, and therefore is optimal. 

The optimality and convergence for the equality problem © can be established similarly to the proof for the 
inequality problem given in lBS02bi l except for the proof for feasibility of the dual variables corresponding to the 
lower and upper bounds (Propositions 8 and 9 in I BS02bl ). Additionally we do not have to prove that p* > 0 
(Proposition 4 in IIBS02bl ) for the equality problem. Therefore we give a complementary proof for the feasibility 
of the dual variables corresponding to the lower and upper bounds. For the equality problem © we have that 
gj(xj) = a,jXj and in the proof we will assume that aj > 0 . 

Lemma 1. Consider the equality problem © and algorithm PIR2 in Section 14.3.11 

(a) If V fe > A k then /A > p k . 

(b) If V fe < A k then p k * < p k . 

Proof. The proof is similar to that of Proposition 7 in llBS02bl . □ 


Proposition 6 (feasibility of dual variables corresponding to bounds). For problem ©, the solution (fl3l > generated 
by PIR2 in Section [43© satisfies the feasibility of the dual variables (p :l and A j) corresponding to the lower and 
upper bounds. 


Proof, (i) For j £ J k U U, we have from (1 1 3cb that pj = 0. 

(ii) For j £ L, we have from (1 1 3cb that pj = f'j(lj) + p k a 3 . We know that all variables Xj with j £ L were 
pegged in the iterations k where V fe > A k . For these iterations k we have x k < lj. Further, from the convexity of 
fj and the assumption that a 3 > 0, we have that p(xj) = —fj(xj)/aj is decreasing in x 3 . Hence, 


Pi 

a,- 




V* = 


~P 


p k * > 0 , 


(14) 


Jjj IU.J Ujj 

where the last inequality follows from Lemma©a). 

(iii) For j £ J k U L, we have from ( 1 1 3cb that A j = 0. 

(iv) For j £ U, we have from ( 1 1 3eb that A j = —<p'Auj) — p k 


* j “j- 


We know that all variables x, with 


j £ U were pegged in the iterations where V ' < A . For these iterations k we have x k > Uj. Further, from the 
convexity of fj and the assumption that a 3 > 0, we have that p(xj) = —f'Axj)/aj is decreasing in x 3 . Hence, 


9i(uj) 


P A 


> - 


g’M) 


k* 


= P k - g k * > o, 


where the last inequality follows from Lemma©b). 


(15) 

□ 


The algorithms in Sections [4. 3. 2114.331 and l4.3.4l will also converge to the optimal solution since they are equiva¬ 
lent to PIR. 


Remark 8 . If we introduce the additional assumption that fj and gj are twice differentiable and that p' > 0, an 
alternative convergence result for the dual relaxation algorithm is found in ISteOll . 
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4.5 Time complexity 


Consider algorithm PIR2 in Section 14.3.11 In Step 0 we need at most 2 n comparisons to determine the primal 
variables and Cqu operations, for a constant Co, to determine if the solution is feasible or not. In Step 1, we solve 
the relaxed problem, which gives C\n operations for a constant C\. In Step 2 we perform at most 2 n comparisons 
to determine the lower and upper sets and we compute A fc and V fe , which gives at most C 2 n operations for a 
constant C 2 . Steps 3.1 and 3.2 give at most 2n + 1 operations. Further, in the worst case the algorithm only pegs 
one primal variable in each iteration, which results in n iterations. Hence, we conclude that the algorithm has a 
complexity of 0(n 2 ). 


5 A quasi-Newton algorithm (NZ) 


Nielsen and Zenios HNZ921 Section 1.4] develop a quasi-Newton method for finding the dual optimal solution //* 
of the problem ©. It is assumed that the objective term (j)j is strictly convex with a derivative </>'■ whose range 
is R. They compare their numerical method with three linesearch methods I HKL80liCL811lTse90l . Their results 
show that their numerical method always performs well compared with the other algorithms. They implement their 
algorithms on a massively parallel computer. This is not our intention. But since the algorithm seems to perform 
well on parallel computers it makes sense to evaluate it on non-parallel computers. 

Let fj, j £ J, be the inverse of </>' such that 

for j £ J, <t>j(fj(fj)) = Ab M G R. (16) 


Similar to ([5]) we conclude that 

Xj(i u) = max{/j, min {fj(a,jn), Uj}}. (17) 

The heart of the algorithm is, like in the breakpoint algorithm, to find /1 such that the primal constraint ( I I b[ i is 
fulfilled. In other words find // such that 



T-( m ) :=b-Y,ajXj 

(m) = 0 . 


(18) 


j&J 




Nielsen and Zenios [NZ92, Section 1.4] define two functions $ 

+ and <l> j : 



^ + (M) := | 

\ min {fj( aj n),Uj}, 

{ ma x{lj,fj( a jn)}, 

O O 

A V 

j G J, 

(19) 

and 





*7M : =| 

I min {fj(ajn),Uj}, 

{ ma x{lj,fj(ajii)}, 

if dj < 0, 
if dj > 0. 

j G J, 

(20) 


Note that for j £ J , if dj > 0 and fj is concave and increasing then <\>f is concave and if a :) > 0 and fj is convex 
and decreasing then <1^ is convex. Further, we define two sets of indices such that J + := { j £ J \ a,j > 0 } and 
J~ := { j £ J | a,j < 0 }. Define two approximation of 4' such that 


jeJ 

= 6— dj min {g (a jff),Uj} — max {lj, g(ajiJ,)}, (21a) 

jeJ+ jeJ- 
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and 


^ 0 ) (//) 

jeJ 

= b ~ Y a :i min {5(ajM), %} - max {(.,■, g(ajfi)}. (22a) 

j'eJ- ieJ+ 

Note that if a ? > 0 and fj is concave then 'f ,+ is convex and if a :l > 0 and fj is convex then is concave. 
Define the sub- and superdifferentials of ’J ,+ , respectively, 'I' - , as 

d'f + (n)~{d£R\{'Z + (n , )--'I ,+ (n)>d(lx'-fi) Vel}, (23a) 

chI/ _ (/z) := { d £ R | — ^“(/z) < d(/jf — /z) V/i'el}. (23b) 

Further, define /z* and x* as the approximate dual and primal solution such that | T (//*) < e where e > 0. The 
algorithm (NZ) follows ( |NZ92l Linesearch 4]): 

Initialization: Set e > 0, k = 0, /z° £ R. 

Iterative algorithm: 

Step 1 (compute step size): 

If T'(/i fc ) > e then 

A fi k+1 := — - where d k £ 9T' + (/x fc ); go to Step 2, 

else if < — e then 

A fi k+1 := — ^ where d k £ 9T'~(/x fe ); go to Step 2, 

else 

Set n* := fi k and determine x* from (IT71) . Stop. 

Step 2 (dual variable update): 

set n k+1 := /. L k + A/r fe+1 . 

Step 3: Set k := k + 1 and go to Step 1. 

The algorithm converges to a value /j,*, such that |\D(//*)| < e if the objective function components 4>j is such that 
the corresponding function T ,+ (/i) is convex or if the corresponding function T “(//) is concave [ NZ92 ] Proposition 
8]. For some problems, the inverse of the derivative might however result in imaginary values. One solution to this 
problem is to consider the equivalent maximization problem of ([2}. he., to maximize^ 

Remark 9. In practice we will choose e > 0 and the algorithm will in most cases stop such that |T'(/r fc )| > 0. 
Hence we will end up with an approximate solution of the optimal dual variable fi*. The map from the dual space to 
the primal might not be linear. Hence the primal error might be larger than we expect, i.e., x* — x* \ » /;* — /r . 
However, there are methods for generating primal optimal solutions from any Lagrangian dual vector (see for 
example I LMOPQ81 ). Another plausible method to find the optimal solution from the approximate solution /i* is 
to use a breakpoint or a relaxation algorithm, starting from /z*. 

6 Method of algorithm evaluation 

This section serves to provide an overview of the procedure for the numerical study. In Section 16.11 we define 
problem instances for the numerical study. Some theory on how the problem instances can be designed follows 
in Section 16721 In Section [6731 we give a brief overview of performance profiles ([ DM021 ) which are used for the 
evaluation of the numerical study. Finally, in Section l674l we describe the computational environment. 
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6.1 Problem set 


For the numerical study we consider five common special cases of problem 0; it also covers the inequality 
problem <03 when //* > 0. Only finite values of the lower an upper bounds d2cl > are considered, i.e., for j £ J, 
lj > —oo and Uj < oo. The five problem cases are briefly specified next: 

Quadratic problem The convex separable quadratic problem is the special case of (0, where 

<f>j{xj ) = - cjXj for j £ J, (24) 

where uij,Cj > 0, j £ J. Numerical studies of algorithms for problem (1241 ) are widely explored; e.g., see 
[ NZ921 Kiw()7 . KiwOS al lKiw08bj . In our numerical study the parameters are randomized such that a,j £ [1, 30], 
Wj £ [1, 20], Cj £ [1, 25], lj £ [0, 3] and Uj £ (3,11], 


Stratified sampling If we have a large population and would like to perform a statistical research among the 
population, it is practically infeasible to examine every single individual in the population. Instead we can stratify 
the population into n strata. An example of a stratum might be people of a certain age. Let M be the number 
of individuals in the entire population and Mj the number of individuals in strata j. If we want to minimize the 
variance of the entire population we need to allocate the number of samples Xj from each strata from: 

<t>o ( x j ) = "j ( y/ for 3 e J ’ (25) 

where uij = Mj /M and pj is an estimate of the variance for strata j. From b in the resource constraint we specify 
the total sample size. In our numerical study the parameters are randomized such that <ij £ [1, 30], rrij £ [5, 30], 
Cj £ [1,4], lj £ [1, 3] and Uj £ (3,15]. 

Sampling According to IIBRS99I , often the objective terms for sampling problems can be written as 

4>j(xj) = Cj/xj for j £ J. (26) 

In our numerical study the parameters are randomized such that a,j £ [1,4], Cj £ [5,30], lj £ [0,3] and u :j £ (3,6]. 

The theory of search In the theory of search problem it is determined how a resource b of time should be spent 
to find an object among n subdivisions of an area with the largest probability. It is assumed that we know the 
probability rrij for an object to be found in subdivision j. The objective component <t>j describes the probability of 
finding the object in subdivision j and takes the form: 

4>j{ x j) = rrij(e~ bjXj — 1) for j £ J. (27) 

This problem is widely explored by Koopman IKoo53l !Koo99fl and the problem is possible to apply to a large 
variation of search problems, e.g., searching for refugees fleeing from Cuba l Sto81 i. In our numerical study the 
parameters are randomized in the following intervals: rrij £ [0.5, —8], bj £ [0.1, 3], cij £ [1,3], lj £ [0, 0.1] and 
Uj £ (0.1, 5]. 

Negative entropy function The negative entropy function mentioned in I NZ92 I: 

<t>j(xj) = Xj log (— - 1) for j £ J. (28) 

cij 

In our numerical study the parameters are randomized such that Cj £ [50, 250], lj £ [20,100] and Uj £ (30, 210]. 
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6.2 Design of problem instances 


Similarly to the numerical study in [ KL98I , we divide our set of test problem instances into groups containing 
different portions of the activities within the lower and upper bounds at the optimal solution. Let H := {j £ 
J | lj < x* < Uj }. Then, the percentage is determined as \H\/n. To motivate this approach, we refer to the 
variance of CPU times for different portions of active activities in lKL98ft . Further, similar to the numerical study 
in lKiw07l . we consider different problem sizes n since the theoretical CPU time for different algorithms vary 
between 0(n) and 0(n 2 ). 

For the problem set, we pseudo randomize the parameters lj, Uj and all the parameters associated with <pj 
and <i :l . In the numerical study we use a linear resource constraint such that gj(xj) = a,jXj, where > 0 for all 
j £ J. This simplifies the design of the problem set, since d> 3 is convex and lj < Uj for all j £ J. We have that 

x*, if ix* = ~$j {x'j)/aj, 

lj, if n* > -<t>'j(lj)/a,j > -(/>'j(uj)/aj, (29) 

Uj, if n* < -<f/j(uj)/aj) < -<t>'j{lj)/aj. 

By using the properties of (l29b we can determine H such that \H\/n = y for any y £ [0,1]. 



6.3 Performance profiles 


Dolan and More BDM02I propose a performance profile for the evaluation of optimization software. The method 
is briefly summarized as follows: Assume that we have a set A of algorithms that consist of n a algorithms and a 
problem set P that consists of n p problem instances. Let t p<a denote the time it takes for algorithm a G A to solve 
problem p £ P. A performance ratio r pa is introduced: 


G*,a(fp,a) - — 


min{f Pi ; | l G A} 


(30) 


For a problem p, the performance ratio is a measure of how fast algorithm a is relative to the fastest algorithm 
solving problem p. Fix a constant t'm such that tm > r>,a for all p G P, a G A and let r p a = xm if algorithm a 
fails to solve problem p. Further, we introduce the distribution 

Pa{r) = — \{p G P I r P} a < t}|, (31) 

n p 


for each algorithm, where | • | denote the cardinality of a set and r G [1, r/y,/]. For algorithm a, the distribution p a 
describes the percentage of problem instances that are solved at least as fast as r times the fastest algorithm for 
problem p. Note that p a { 1) is the percentage for algorithm a being the fastest. Moreover, lim,.^,^ p s (r) is the 
probability that algorithm a will solve a problem p in P. If we have a large problem set P then p„(r) will not be 
affected much by a small change in P ( IDM021 Theorem 1]). 


6.4 Program language, computer and code 

The algorithms are implemented in Fortran 95, compiled with gfortran under mac OS X 10.8.2 (2.5 GHz Intel 
Core i5, 4 GB 1600 MHz DDR3). 


7 Computational experiments 

In this section we present the results from numerical experiments of the problems defined in Section [6TI If nothing 
else is mentioned, 100 problem instances of each problem is evaluated for each problem size. Further, the problem 
instances are designed such that different values of \H\/n are considered, see Section [(PI 
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The development of numerical experiments for the problem (1) is illustrated in Table [T] where we cite the size 
of the largest test problem reported during each decade; the table is an extension of Table 2 in |PatQ81l . 


TABLE 1: Largest problem instances solved for each algorithm class through the decades. 


Decade 

50s 60s 

70s 

80s 

90s 

00s 


10s 

Breakpoint algorithm 

2 60 

12 

200 

10 4 

2-10 6 

30 

10 6 

Relaxation algorithm 

- 

200 

200 

10 4 

2-10 6 

110 

10 6 


In Section [7711 we present performance profiles for the relaxation algorithms defined in Section 0] In Section 
17.21 we evaluate the pegging process for both the breakpoint and the relaxation algorithm defined in Section 0 
In Section POl the best performing relaxation algorithm and breakpoint algorithm are compared with the quasi- 
Newton algorithm in Section [5] In Section 17741 we give a critical review of the numerical study. Finally, in Section 
[8j we make some overall conclusions. 

In the following Figures [2^4] we show performance profiles, as defined in Section 16.31 for three competing 
sets of algorithms. To summarize the appearance of these plots, for each plot and corresponding competing set 
of algorithms, the graph for one algorithm shows the portion (between 0 and 1 on the y-axis) of the problem 
instances considered that are solved within r (on the x-axis) times the fastest algorithm in the test for problem p. In 
particular, the value of p a (l) is the portion of the problems in which algorithm a is the fastest, and lim T _» rM p a {j) 
is the probability that algorithm a will solve a problem p in P. (The latter information is particularly illustrative 
and relevant for Figure 0]) 

7.1 Evaluation of the relaxation algorithm 

We evaluate the different paths in Figure Q] First, recall that PIR2 determines the primal variables while DIR2 
determines the dual variable. As can be seen in the leftmost performance profile in Figure[2] DIR2 is the fastest in 
67.8% of the problem instances solved. The result is what we can expect from theory, since DIR2 does not need 
as many operations in Step 1; see Sections [4. 3. II and 14. 3. 21 PIR2 is faster in 32.2% of the problems solved; the 
latter cases stem mainly from to the negative entropy problem (1281) but also from the theory of search problem (f27l) 
and the stratified sampling problem (l25l > when n and \H\/n is small. The reason for PIR2 to be faster in almost 
all cases for the negative entropy problem is due to the computationally simple expression of the primal variables 
x k (p k ) = ^ -while the dual variable is evaluated from u k = logV. ^ Tk — log 5. Also the computations 

J 2^/j£jk c j *■— J 

of the breakpoints might be a larger part of the total time when the equations for the primal/dual variables are 
simple. 

We recall that DER2 applies explicit evaluation, DIR2 applies implicit evaluation and DBR2 applies the the¬ 
oretically most profitable evaluation. The results are in favour of DBR2 and DIR2; see the rightmost performance 
profile in Figure[2] DBR2 is fastest in 70.6%, DIR2 is fastest in 25.9% and DER2 is fastest in 4.7% of the problems 
solved. (Note that 70.6% + 25.9% 4- 4.7% = 100.2%; in some cases two algorithms are equally fast.) 

Figure[2]also shows that the performance of DIR2 and DBR2 are quite similar. For just a few cases DIR2 is 
more than 1.30 times slower than the fastest algorithm while for as few cases DBR2 is only 1.2 times slower than 
the fastest algorithm. On average DBR2 performs somwhat better than DIR2. 

Conclusion For the problem set considered, it is more profitable to evaluate the dual variable even if we have to 
compute all the breakpoints at the beginning of the algorithm. Hence for the problem set considered DIR2 seems 
to outperform PIR2. 
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FIGURE 2: The leftmost figure shows the performance profiles for PIR2 and DIR2 and the rightmost for DER2, DIR2 and 
DBR2. The numerical experiment was done according to the description in Section [6~2l (or n = 50 • 10 3 ; 100 ■ 
10 3 ;200 • 10 3 ;500 • 10 3 ; 10 6 ; 2 • 10 6 and for the problems in Section [6~i1 . The algorithms were implemented 
uniformly. 


For the problem set considered, the results show that in most cases it is more profitable to evaluate a solution 
x fe implicitly. The small difference between the performance of DBR2 and DIR2 implies that in most cases 
| J k \ > 2|U U {^ k )|. We conclude that DBR2 performs slightly better than DIR2 which agrees with the 
theory in Section !?. 3. 41 


7.2 Evaluation of the pegging process 

We compare the three pegging approaches for the breakpoint and relaxation algorithms described in Sections [3.31 - 
13.51 and 14.3.51 respectively. For the breakpoint algorithm median search is implemented similarly to I FTVF92 
Section 8.5]. 

The performance profiles in FigureOshow that 5-sets pegging is nearly always the fastest (MB5 in 95% and 
DBR5 in 94% of the problem instances solved). Additionally, when 5-sets pegging is not fastest it is never more 
than 10% slower than the fastest algorithm. It is obvious that for the few cases when 2-sets pegging performs best 
are when few optimal primal variables equals the lower or upper bounds (\H\/n is small). This follows from the 
fact that 2-sets pegging does not check if the variables belong to the additional sets that are used in 3- and 5-sets 
pegging. 


Conclusion In general 5-sets pegging is the most profitable. 


Remark 10. Numerical experiments for a breakpoint algorithm applying bisectional search of a sorted sequence 
of breakpoints, using quicksort as in | Knu98 , Section 5.2.2], was also performed. The algorithm was on average 
1.5 times slower than MB2. Also the algorithm in | Zip80| , where the upper bounds are relaxed, was evaluated. 
The algorithm performed poorly, on average 2.9 times as slow as the breakpoint algorithm applying sorting. 


7.3 A comparison between relaxation, breakpoint, and quasi-Newton methods 

In this section we compare the best performing relaxation algorithm (DBR5), the best performing breakpoint al¬ 
gorithm (MB5), and the numerical quasi-Newton method (NZ) in Section^ For NZ we use the stopping criterion 
a 3 x i _ j | < o.Ol which is weaker than the stopping criterion in 1 NZ921 (< 10~ 4 ). This will of course give us 
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FIGURE 3: The leftmost figure shows the performance profile for MB2, MB3 and MB5 and the rightmost for DBR2, DBR3 
and DBR5. The numerical experiment was done according to the description in Section [6!2l for n = 50 • 10 3 ; 100 • 
10 3 ; 200 ■ 10 3 ; 500 ■ 10 3 ; 10 6 ; 2 • 10® and for the problems in Section ETT1 The algorithms were implemented 
uniformly. 

an approximate optimal solution only. The initial value of the dual variable // is set to the mean of the breakpoints 
fi° = (Ij)/ a j + ( t > j( u j)/ a j)/(^ n )- ^ If 16 algorithm does not converge within 100 CPU seconds, 

we start over with the mean of the breakpoints corresponding to the lower breakpoints /r° = ‘A? (h)/ a j)/ n - 

Similarly, if the algorithm does not converge within 100 CPU seconds then we start over with the mean of the 
breakpoints corresponding to the upper breakpoints: /i° = / a j) ! n - If the algorithm does not termi¬ 

nate within 300 CPU seconds then we terminate the algorithm and consider the problem as unsolved. 

The reader may have noticed that the breakpoint algorithms in Sections [3.3113. 5l are given without the use of 
the sets for the pegging process ( L , U, M, L_, [/+) while the relaxation algorithms in Section l4.3.1H4.3.5l are given 
with the use of these sets. Of course the notation using these sets, as for the relaxation algorithms, are theoretically 
more elegant. However, our implementations of the algorithms in Fortran 95 showed that for breakpoint algorithms 
it was in general more profitable to not use the pegging sets. 

Figure 0] shows the performance profile for the algorithms. In general, we can see a significant advantage of 
using DBR5 since it is fastest for 99.0% of the problem instances solved. Moreover, when DBR5 is not the fastest 
algorithm it is almost as fast as the fastest. From the figure we also see that MB5 is more than 2.7 times slower 
than the fastest algorithm in 50% of the problem instances solved. 

The Newton method terminated the fastest for only 1% of the test problems, which probably is due to good 
initial values. Also we should have in mind that NZ terminates with an approximate solution only. The Newton 
method performs relatively well for the quadratic problem (l24l> . the theory of search problem ( |27} . and the negative 
entropy problem ( 128} when \H\/n > 0.8. It is not very successful for the stratified sampling problem ([25} and the 
sampling problem (126} . We can see from the performance profile that in 28% of the problem instances tested, NZ 
is more than 5.5 times slower than DBR5 (the fastest algorithm). In 5.3% (158/3000) of the problem instances the 
algorithm does not solve the problem; these problem instances mainly stem from the stratified sampling problem 
(125} and the sampling problem (126} when \H\/n < 0.3. 

By comparing the performance profiles in Figure 0] we can see that they are very similar. In other words, 
relative to each other, the performance of the algorithms seems to be similar for different problem sizes n. 

Considering the breakpoint algorithm MB5, Figure 0] shows that the CPU time increases when the number of 
optimal primal variables strictly within the bounds \H\/n grows for the problems in the problem set. Concerning 
the relaxation algorithm DBR5 the result shows no significant dependence of \H\/n, neither does NZ. 
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FIGURE 4: Performance profile for DBR5, MB5 and NZ. (Left) The numerical experiment was done according to the descrip¬ 
tion in Section l6^2l for n = 50 ■ 10 3 , 100 • 10 3 , 200 ■ 10 3 , 500 ■ 10 3 , 10 6 , and 2 ■ 10 6 for the problems in Section l6Tl 
(Right) The numerical experiment was done according to the description in Section [6^21 but with 10 problem in¬ 
stances and for n = 4 ■ 10 6 , 6 ■ 10®, 8 ■ 10®, 10 • 10®, 15 • 10®, 20 ■ 10®, 25 • 10®, and 30 • 10® for the problems in 
Section l6.ll 



Average |H|/n 


FIGURE 5: The average CPU time for n = 10® is plotted as a function of the average portion of the optimal variables that 
obtain a value between the lower and upper bounds. 

In Figure[6] we can see that the CPU time for MB5 and DBR5 is linear for n G [50 • 10 3 , 2 • 10 6 ] respectively 
n € [50 • 10 3 ,30 • 10 6 ]. This is impressive, since DBR5 has a worst-case time complexity of 0(n 2 ). For problem 
sizes a bit larger than 30 • 10 6 variables the CPU time tends to increase faster. However, additional computations 
on a more powerful computer have shown that, for the relaxation algorithm, the quadratic problem has a linear 
time complexity up to n = 90 ■ 10 6 , the sampling problem up to n = 100 ■ 10 6 , the theory of search problem up to 
n = 70 • 10 6 and the negative entropy problem up to n = 110 • 10 6 . Hence we would also like to stress that the 
linearity probably is constrained by the memory of the computer rather than the algebra in the algorithm. 

Tables [2^6] show the mean CPU time in seconds, for solving the problems in the problem set for DBR5, NZ 
and MB5. The tables show the great advantage of using DBR5. 
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FIGURE 6: The average CPU time is plotted as a function of the size of the problem n for DBR5 and MB5. The values in 
the left figure are the mean of 100 problem instances solved and the values in the right figure are the mean of 10 
problem instances solved. 

Table 2: The average CPU times algorithms MB5, DBR5 and NZ for solving the quadratic problem d24l l. Each value is the 
mean of 100 randomized computations and the CPU times are given in seconds. 


n 

50,000 

100,000 

200,000 

500,000 

1,000,000 

2,000,000 

DBR5 

0.0031 

0.0065 

0.0136 

0.0346 

0.0702 

0.1438 

NZ 

0.0083 

0.0179 

0.0378 

0.1018 

0.1923 

0.3822 

MB5 

0.0097 

0.0205 

0.0428 

0.1144 

0.2346 

0.4807 


Table 3: The average CPU times algorithms MB5, DBR5 and NZ for solving the stratified sampling problem ( 125b . Each value 
is the mean of 100 randomized computations and the CPU times are given in seconds. For NZ only the cases when 
NZ found the optimal solution have been considered. 


n 

50,000 

100,000 

200,000 

500,000 

1,000,000 

2,000,000 

DBR5 

0.0040 

0.0082 

0.0167 

0.0436 

0.0875 

0.1773 

NZ 

0.0379 

0.0952 

0.1529 

0.3560 

0.9026 

1.8299 

MB5 

0.0098 

0.0214 

0.0438 

0.1178 

0.2390 

0.4981 


Table 4: The average CPU times algorithms MB, DBR5 and NZ for solving the sampling problem ( 126b . Each value is the 
mean of 100 randomized computations and the CPU times are given in seconds. For NZ only the cases when NZ 
found the optimal solution have been considered. 


n 

50,000 

100,000 

200,000 

500,000 

1,000,000 

2,000,000 

DBR5 

0.0032 

0.0069 

0.0146 

0.0375 

0.0788 

0.1597 

NZ 

0.1069 

0.2403 

0.5670 

1.4319 

3.8036 

6.3949 

MB5 

0.0089 

0.0192 

0.0397 

0.1066 

0.2213 

0.4522 
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Table 5: The average CPU times algorithms MB5, DBR5 and NZ for solving the theory of search problem ( 127b . Each value 
is the mean of 100 randomized computations and the CPU times are given in seconds. 


n 

50,000 

100,000 

200,000 

500,000 

1,000,000 

2,000,000 

DBR5 

0.0063 

0.0132 

0.0259 

0.0679 

0.1382 

0.2778 

NZ 

0.0187 

0.0434 

0.0940 

0.2669 

0.4511 

0.9258 

MB5 

0.0127 

0.0273 

0.0549 

0.1454 

0.2983 

0.6090 


Table 6: The average CPU times algorithms MB5, DBR5 and NZ for solving the negative entropy problem ( 128b . Each value 
is the mean of 100 randomized computations and the CPU times are given in seconds. 


n 

50,000 

100,000 

200,000 

500,000 

1,000,000 

2,000,000 

DBR5 

0.0036 

0.0074 

0.0157 

0.0406 

0.0837 

0.1698 

NZ 

0.0071 

0.0171 

0.0310 

0.0819 

0.1709 

0.3191 

MB5 

0.0096 

0.0206 

0.0427 

0.1131 

0.2326 

0.4774 


Summary Even if NZ only finds an approximate solution and even if the stop criterion is weak, NZ is outper¬ 
formed by DBR5. Hence, we see no reason to evaluate methods for finding the exact optimal solution x* from the 
approximate solution x* generated by NZ. In many cases, NZ fails to solve the stratified sampling problem (j25j 
and the sampling problem (l26l >. This is probably due to the stiffness of the problems. 

MB5 performs very well for problems where \H\/n is small. This is due to the pegging process: since we 
peg a large part of x k € J k in each iteration, the problem is reduced in the next iteration. This follows from the 
reduction of breakpoints by half in each iteration; note that this does not hold for the relaxation algorithm. This 
might be the reason for MB5:s dependency of | Il\/n. However, DBR5 performs best for almost all of the problem 
instances solved and it is worth noting that DBR5 never performs poorly. 

7.4 A critical review 

When we set up the initial values for NZ, we tried different initial values and concluded that the one used was the 
most profitable on average. However, especially if we customize initial values for each problem, there may still be 
potential improvements available. 

Considering the results in Section 177X1 we consider the dual algorithm (DIR2) to outperform the primal algo¬ 
rithm (PIR2) since it performs better in general. Then we continued to evaluate several modification of the dual 
algorithm while similar modifications of the primal relaxation algorithm were not evaluated, i.e., we don’t evaluate 
blended evaluation or 3- and 5-sets pegging for the primal algorithm. However since the plausible modification of 
the primal algorithm would take the same form as for the modification of the dual algorithm it is assumed that a 
modified dual algorithms will outperform a similar modified primal algorithm. 

8 Conclusion 

We have complemented the survey in HPat08l on the resouce allocation problem at hand, and introduced, and 
critically evaluated, new implementations of breakpoint and relaxation algorithms for its solution. 

The results show that our new implementations (DIR2 and DBR2) of the relaxation algorithm outperform 
the earlier algorithms (PIR2 and DER2). Hence we should evaluate the dual variable for the relaxation algorithm. 
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i.e., DIR2 outperforms PIR2. Moreover, it is more profitable in theory, as well as in practice, to apply blended 
evaluation, i.e., DBR2 outperforms DER2 and DIR2. Our results, as well as the results in [ RJL92. KL981lKTw08bl . 
imply that the relaxation algorithm is to prefer when a closed form of the dual variable /i can be found. 

We introduced 3- and 5-sets pegging for the relaxation algorithm and showed that it is most profitable to apply 
5-sets pegging, which also holds for the breakpoint algorithm. 

For the problems considered, MB5 and DBR5 have a practical time complexity of 0(n) also for very large 
values of n. We also showed that the relaxation algorithm DBR5 performs better than both the Newton-like 
algorithm NZ and the breakpoint algorithm MB5. Potential future improvements include the implementation of a 
pegging method for a Newton-like algorithm and/or a hybrid of the different algorithms (NZ, DBR5 and MB5), 
and, hence, it would be of interest to compare the best algorithms from our study to these. 

The findings made herein can most certainly be profitably utilized also in the efficient solution of the more 
complex versions of the resource allocation problem discussed, for example, in the books | Mje83[llK88 , Lusl2; l. 
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